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ABSTRACT 


This  thesis  investigates  and  compares  time  and  wavelet- domain  denoising  tech¬ 
niques  where  received  signals  contain  broadband  noise.  We  consider  how  time  and 
wavelet- domain  denoising  schemes  and  their  combinations  compare  in  the  mean  squared 
error  sense.  This  work  applies  Wiener  prediction  and  Median  filtering  as  they  do  not 
require  any  prior  signal  knowledge.  In  the  wavelet- domain  we  use  soft  or  hard  threshold 
on  the  detail  coefficients.  In  addition,  we  explore  the  effect  of  these  wavelet- domain 
thresholding  techniques  on  the  coefficients  associated  with  cycle- spinning  and  the  newly 
proposed  recursive  cycle -spinning  scheme.  Finally,  we  note  that  thresholding  does  not 
make  an  attempt  to  de- noise  coefficients  that  remain  after  thresholding;  therefore  we 
apply  time  domain  techniques  to  the  remaining  detail  coefficients  from  the  first  level  of 
decomposition  in  an  attempt  to  de- noise  them  further  prior  to  reconstmction.  This  thesis 
applies  and  compares  these  techniques  using  a  mean  squared  error  criterion  to  identify 
the  best  performing  in  a  robust  test  signal  environment.  We  find  that  soft  thresholding 
with  Stein’s  Unbiased  Risk  Estimate  (SURE)  thresholding  produces  the  best  mean 
squared  error  results  in  each  test  case  and  that  the  addition  of  Wiener  prediction  to  the 
first  level  of  decomposition  coefficients  leads  to  a  shghtiy  enhanced  performance. 

EinaUy,  we  illustrate  the  effects  of  denoising  algorithms  on  longer  data  segments. 
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XIV 


EXECUTIVE  SUMMARY 


This  thesis  investigates  time  and  wavelet- domain  denoising  algorithms  in  a  robust 
signal  environment,  and  performs  a  comparative  analysis  between  the  various  denoising 
schemes.  Wavelet-based  denoising  has  mihtary  apphcations  in  acoustic  signal  process¬ 
ing,  surveillance  and  reconnaissance,  and  various  communication  apphcations.  This 
work  reviews  the  basic  theory  of  wavelet  denoising  leading  to  the  discrete  orthogonal 
wavelet  transform.  In  addition,  we  present  and  apply  cycle -spinning  denoising,  which 
serves  to  remove  the  negative  effects  of  translation  variance  found  in  the  orthogonal 
wavelet  transform.  FinaUy,  we  investigate  recursive  cycle- spinning  and  a  combined 
time- wavelet  denoising  technique  using  a  windowed  Wiener  predictor. 

Results  show  that  soft  thresholding  of  the  translation- invariant  wavelet  transform 
coefficients  with  the  Stein’s  Unbiased  Risk  Estimate  (SURE)  threshold  outperforms  the 
other  schemes  considered  in  Mean  Squared  Error  (MSE)  performance  for  the  shorter 
signal  lengths  for  the  signals  investigated.  Eor  the  longer  data  length  experimental 
signals  considered  the  translation- invariant  visual  threshold  showed  the  best  denoising 
abihty.  In  addition,  no  difference  could  be  seen  in  the  performance  between  orthogonal 
wavelet  and  translation- invariant  schemes  for  the  long  data  sets  considered.  Results 
further  show  recursive  cycle -spinning  using  the  SURE  threshold  for  smaller  data  sets  and 
the  visual  threshold  for  longer  data  sets,  though  not  yet  optimized,  have  excellent  poten¬ 
tial  as  a  denoising  alternative.  In  addition,  we  found  through  experimentation  that  the 
window  size  of  the  Wiener  predictor  requires  a  priori  information  for  optimal  perform¬ 
ance. 
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1.  INTRODUCTION 


A.  THESIS  OBJECTIVE 

This  study  considers  a  frequency- based  approach  to  denoising.  The  noisy  signal 
is  transformed  into  the  frequency  domain  using  wavelet  transforms  where  denoising  is 
conducted.  Next  the  denoised  signal  is  obtained  by  back- transforming  the  processed 
information  to  the  time  domain.  Further  developments  in  this  area  will  lead  to  advances 
in  range  and  sensitivity  with  regard  to  acoustic  signal  processing,  surveillance  and  re¬ 
connaissance,  and  various  communications  apphcations. 

Donaho  and  Johnstone  [1,2]  introduced  several  methods  where  information  in  the 
wavelet  transform  domain  may  be  used  to  de- emphasize  the  noise- contaminated  fre¬ 
quency  contributions  prior  to  performing  the  inverse  transform,  thereby  producing  a 
cleaner  system  output.  This  thesis  explores  these  methods  and  compares  resulting 
performances.  The  goal  of  this  thesis  is  to  determine  which  decomposition  type  and 
combination  of  signal  processing  schemes  achieves  the  cleanest  signal  output. 

Figure  1.1  presents  the  schematic  overview  of  the  concepts  considered  in  this 
work.  It  methodically  outiines  and  organizes  this  work  into  six  chapters  from  the  theory 
behind  the  work  to  the  results  and  conclusions.  Chapter  II  introduces  the  wavelet  trans¬ 
form  and  illustrates  its  advantages  over  Fourier  analysis.  Chapter  IH  introduces  signal 
and  noise  elements  and  the  associated  assumptions  that  went  into  producing  robust  test 
signals.  Chapter  IV  first  presents  various  wavelet  thresholding  techniques.  Next,  it 
introduces  translation  invariance  concepts  and  shows  how  they  are  obtained  via  processes 
called  cycle- spinning  and  recursive  cycle- spinning.  Finally,  it  introduces  a  windowed 
Wiener  prediction  technique  and  includes  a  low  computational  cost  scheme  called 
median  filtering.  Chapter  V  applies  these  techniques  to  test  data.  Chapter  VI  presents 
conclusions  and  recommendations. 
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B.  THESIS  ORGANIZATION 


Figure  1.1:  Thesis  outline  flow  structure. 
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11.  MULTIRESOLUTION  ANALYSIS 


A.  INTRODUCTION 

This  chapter  presents  an  overview  of  Multiresolution  Analysis  (MRA)  as  it  leads 
to  an  important  family  of  signal  processing  tools.  The  signal  processing  apphcations 
rooted  in  MRA  have  the  distinct  abihty  to  span  both  time  and  frequency,  as  multiresolu¬ 
tion  analysis  may  be  viewed  as  an  extension  of  the  Short  Time  Fourier  Transform 
(STFT). 


I.  Short  Time  Fourier  Transform 

The  Fourier  Transform  is  a  widely  used  transformation  and  is  defined  for  a  con¬ 
tinuous  signal  f(t)  as: 

+00 

F(a))  =  j (2.1) 

The  complex  basis  function,  also  referred  to  as  a  kernel,  in  Equation  (2.1) 
operates  on  /  (t)  over  aU  time,  creating  a  one- dimensional  representation  of  the  signal  in 
the  frequency  domain.  Note  that  it  is  often  desirable  to  locahze  the  specific  time  at  which 
specific  signal  characteristics  occurred;  however,  there  is  no  allotted  time  dimension  in 
the  Fourier  Transform.  Time  dependency  capabihty  is  introduced  in  the  Short  Time 
Fourier  Transform  (STFT)  with  a  predetermined  fixed  duration  shding  window  g(t-x) 
centered  at  t .  The  resulting  two-dimensional  SIFT  of  the  signal  /  (t)  is  defined  as: 

+00 

F((D,x)=  j/(x)g(x-t)e-^‘^Jx.  (2.2) 

The  signal  /  (t)  and  the  Fourier  Transform  F((o)  are  contained  entirely  in  one¬ 
dimensional  space,  whereas  the  STFT  transform  F((0,x )  is  represented  in  a  two-dimen¬ 
sional  space.  The  magnitude  squared  of  the  STFT  is  called  the  Spectrogram  as  shown  in 
Figure  2.1,  which  provides  a  time-localized  frequency  content  of  the  transformed  signal 
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/  (t) .  Figure  2. 1  illustrates  signal  presence  in  red  and  a  lack  of  signal  presence  in  blue, 
obtained  for  the  following  constant  amplitude  linear  chirp: 


/(0  =  sin 


271 

V  ^ 


uo 


(2.3) 


where  is  set  to  8000  Hz,  and  the  frequency  /„(0  varies  linearly  from  1  to  2000  Hz. 


Time  (n) 


Figure  2. 1 :  Linear  chirp  spectrogram. 

The  frequency  resolution  with  a  fixed  duration  sliding  window  g(t-x)  is  fixed. 

In  order  to  get  a  perfect  picture  of  the  signal  space  from  the  Spectrogram,  it  would  be 
desirable  to  have  infinite  time  and  frequency  resolution.  Unfortunately,  the  STFT  time 
and  frequency  resolutions  are  limited  by  the  Heisenberg  Uncertainty  Principle  [3]: 

(2.4) 

4 
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^  2 

where  a ,  and  o  ^  represent  the  temporal  variance  and  frequency  variance,  respectively. 

The  best  case,  with  equahty  in  Equation  (2.4),  is  found  when  using  a  Gaussian  window 
git-x  )  in  the  STEP.  This  results  in  the  Gabor  Transform  [4,5] .  For  discrete  time 
signals  a  shortened  window  leads  to  fewer  samples,  and  with  fewer  samples  the  resultant 
frequency  space  will  have  fewer  bins.  Hence,  intuitively,  frequency  resolution  is  lost  as 
time  becomes  more  locahzed. 

Note  that  time- frequency  characteristics  may  change  with  time  when  a  signal  is 
transmitted  through  various  channels,  or  the  signal  or  channel  is  time- varying.  As  a 
result,  a  small  time  window  is  necessary  and  frequency  resolution  is  degraded  when  the 
frequency  content  of  a  signal  or  the  channel  characteristics  change  rapidly  with  time.  A 
larger  time  window  may  be  sufficient  (resulting  in  better  frequency  resolution)  when  the 
signal  frequency  characteristics  change  slowly.  For  unknown  signal  characteristics  an 
array  of  window  sizes  is  required  in  order  to  find  the  ideal  time  window  and  frequency 
resolution  combination.  Therefore,  the  fixed  window  size  of  the  STFT  limits  its  abihty  to 
span  both  time  and  frequency  of  unknown  signals  with  resolution  well  matched  to  the 
signal  characteristics. 

From  the  above  discussion,  a  transform  method  that  does  not  utihze  a  fixed  win¬ 
dow  appears  attractive,  and  Multiresolution  Analysis  (MRA)  techniques  do  not  suffer 
from  such  limitations.  Next,  we  discuss  MRA  techniques,  which  do  not  have  a  fixed 
window  constraint. 

2.  Multiresolution  Theory 

In  much  the  same  way  as  the  STFT,  MRA  requires  an  operator  to  project  the  sig¬ 
nal  /  (0  into  another  domain  or  vector  space.  One  of  the  advantages  is  that  the  MRA 
operator  does  not  use  a  fixed  window  size. 

The  Venn  Diagram  in  Figure  2.2  is  meant  to  show  that  vector  space  Vj  is  a  lower 
resolution  approximation  of  Fet  Wj  represent  the  loss  of  information  due  to  the 
approximation,  then  the  vector  sum  of  Wj  and  Vj  fuUy  constitutes  Vj.i.  Further,  there  is 
no  overlap  or  repeat  of  information  if  Wj  and  Vj  are  orthogonal.  Consider  the  entire 
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frequency  domain  from  0  to  p  as  the  vector  space  Vj+j,  as  shown  in  Figure  2.3.  Then, 
part  of  the  frequency  domain  is  spanned  by  Wj  and  the  other  part  is  spanned  by  Vj.  A 
natural  extension  of  the  above  MRA  principles  is  shown  as  follows  [6]: 

Vj.i  =  y,©w,©w,,,©--©Wj 

Appox  Detail 


Figure  2.2:  Venn  diagram  illustrating  MRA  basic  principle  of  successive  approxima¬ 

tions.  After  Ref  [7]. 


Figure  2.3:  Frequency  spectmm  divided  into  multi- resolution  bands.  After  Ref  [7]. 

MaUat  [3,5]  found  that  MRA  could  be  implemented  with  minimal  redundancy  by 
employing  orthonormal  filters  consisting  of  a  lowpass  filter  and  multiple  bandpass  filters. 
The  lowpass  filter  provides  the  approximation  coefficients  and  the  bandpass  filters 
provide  the  detail  coefficients. 

Figure  2.4  shows  the  main  difference  between  MRA  and  the  STFT  is  in  the  filter 
bandwidths.  Another  promising  attribute  of  MRA  is  that  many  different  conventional 
filters  can  be  used  to  implement  it.  An  additional  advantage  is  that  the  operator  itself  can 
be  real,  and  thus  the  coefficients  will  be  real  resulting  in  a  real- valued  transform.  Simple 
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(2.5) 


operations  such  as  the  apphcation  of  a  Finite- duration  impulse  response  (FIR)  filter  can 
be  applied  to  obtain  real- valued  transform  coefficients  if  needed. 

The  filters  chosen  to  form  the  detail  and  approximation  coefficients  are  called 
wavelet  and  scaling  functions  respectively.  These  functions  will  be  explored  in  continu¬ 
ous  and  discrete  time  in  the  following  section. 


(a)  (b) 

Figure  2.4:  Frequency  spectrum  partitioning  for:  (a)  MRA;  (b)  S  lFl  . 

B.  WAVELET  TRANSFORM 


Rather  than  developing  a  different  bandpass  filter  for  each  scale,  MRA  theory  al¬ 
lows  for  a  scaled  and  translated  “mother”  basis  function  or  wavelet.  Wavelet  basis 
functions  act  in  much  the  same  way  as  the  kernel  in  the  Fourier  Transform  and  are 
orthogonal,  normalized  to  unity,  have  finite  time  duration  or  compact  support,  and  their 
area  sums  to  zero. 


1.  Continuous  Wavelet  Transform 

The  continuous  wavelet  transform  (CWT)  is  similar  to  the  Fourier  Transform  in¬ 
troduced  in  Equation  (2.1).  The  wavelet  is  the  basis  function  or  operator  for  the  wavelet 
transform,  as  e  was  the  operator  for  the  Fourier  Transform. 

The  equation  for  the  wavelet  basis  function  at  scale  or  resolution  a  and  transla¬ 
tion  b  of  the  mother  wavelet  v|/  (t)  is  defined  as  [3,8,9]: 
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/  N  1  f t-b'\ 

■sja 


V  «  y 


(2.6) 


Figure  2.5  (a)  shows  how  the  resolution  factor  a  is  used  to  normalize  and  scale 
the  mother  wavelet.  Note  that  the  upper,  middle,  and  lower  figures  in  Figure  2.5  (a) 
illustrate  the  effect  on  the  wavelet  when  a  is  0,  1,  and  2,  respectively.  This  shows  how 
the  wavelet  is  “squeezed”  in  time  as  a  is  increased.  Figure  2.5  (b)  demonstrates  the 
effect  of  the  translation  factor  b,  where  the  upper,  middle,  and  lower  figures  demonstrate 
the  effect  on  the  wavelet  when  b  is  0, 1,  and  2,  respectively.  Hence,  the  wavelet  is 
translated  to  the  right  in  time  as  b  is  increased.  Note  above  that  with  the  CWT  the 
resolution  factor  a  can  be  chosen  to  be  any  positive  real  number,  and  the  translation  b 
can  be  any  real  number.  The  CWT  coefficients  are  given  by: 


da,b  =  I  f(ty\fl,b(Odt. 


(2.7) 


Parseval’s  Theorem  states  that  the  energy  contained  in  a  signal  in  the  time  domain 
is  the  same  as  in  another  domain  if  the  coefficients  are  determined  by  applying  a  set  of 
orthonormal  basis  functions  v(t)  [6],  leading  to: 


f  /(t)v(0dt  =  —  ?  (2.8) 

L  27t  _i 

Note  that,  Parseval’s  theorem  applies  both  to  the  orthogonal  wavelet  and  Fourier 
Transform  operations  since  they  both  use  orthogonal  basis  sets.  Thus,  the  CWT  coeffi¬ 
cients  may  be  expressed  in  terms  of  the  frequency  information  using  Parseval’s  relation¬ 
ship,  which  leads  to: 

=  (2.9) 

Parseval’s  relation  helps  demonstrate  that  the  CWT  coefficients  are  representa¬ 
tive  of  signal  energy  and  larger  coefficients  correspond  to  greater  signal  energy.  The 
wavelet  is  useful  in  data  compression  as  well  because  relatively  few  coefficients  hold 


8 


nx)st  of  the  signal  energy  [1,2, 3, 5, 8].  This  attribute  is  also  very  valuable  to  the  art  of 
denoising. 
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Figure  2.5:  Haar  wavelet:  (a)  effect  of  scaling;  (b)  effect  of  translation. 


Recall  that  the  STFT  has  fixed  time  and  frequency  resolution  as  the  time  window 
length  is  fixed,  while  the  CWT  does  not  have  such  constraints.  The  CWT  time  window 
size  gets  continuously  smaller  as  the  scale  increases,  while  the  frequency  resolution  gets 
larger. 

Wavelets  come  in  different  shapes  and  sizes  and  can  be  chosen  to  meet  situational 
requirements.  Figure  2.6  shows  examples  of  popular  multipurpose  wavelet  basis  func¬ 
tions.  Notice  that  wavelets  have  sharper  transient  characteristics  than  the  sinusoidal 
based  Fourier  Transform  basis  function.  Thus,  the  wavelet  is  better  suited  for  detecting 
transient  signals,  whereas  the  STFT  is  ideal  for  sinusoidal  signals. 
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Figure  2.6:  Assortment  of  wavelets.  After  [10]. 

The  CWT  is  theoretically  significant;  however  the  Discrete  Wavelet  Transform 
(DWT)  is  amenable  to  fast  digital  signal  processing  apphcations.  The  DWT  takes 
discrete  time  input  samples  and  discrete  translation  and  scale  factors,  and  eliminates  high 
processing  costs  associated  with  redundancy  found  in  the  CWT.  The  DWT  is  discussed 
further  in  the  following  section. 

2.  Discrete  Wavelet  Transform 

The  DWT  is  a  clever  discretization  of  the  CWT  that  enables  a  digital  implemen¬ 
tation  of  the  CWT.  From  Equation  (2.5)  and  its  accompanied  discussion  the  approxima¬ 
tion  coefficients  and  detail  coefficients  can  be  written  in  terms  of  the  scaling  function 
(j)(n  -k)  and  the  wavelet  function  \|/  {2^ n-k)  respectively,  as  shown  below  in  Equation 
(2.10).  The  wavelet  function  is  the  highpass  filter,  while  the  scaling  function  is  the 
lowpass  filter.  Orthogonahty  between  all  the  basis  functions  at  various  scales  and 
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translations  prevent  redundancy.  The  fundamental  DWT  expansion  of  the  discrete  signal 
f{n),  where  a  =2^  is  given  by: 


f  in)  =  Co  (« -  ^ )  +  E E  2'' V  [l^n-k). 

k  k  7=0 


(2.10) 


Utihzing  MRA  principles,  Mallat  [3]  proposed  that  within  the  scope  of  MRA  a 
filter- bank  of  highpass  and  lowpass  filters,  gin)  and  hin),  to  make  use  of  coefficients  of 
higher  resolution  to  synthesize  the  lower  resolution  coefficients  as  shown: 


'  J+l,k 


''  n  3 

- k 

V  ^  J 


L* 


dj+i,k 


^  n  3 

- k 

V  ^  J 
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'j,k- 


(2.11) 


The  following  reconstmction  equation  follows  as: 


^j,k  =  E  h{k-2n)cj^^  ,+  g{k-2n)dj^,„  (2.12) 

«=-oo  «=-oo 

where  hin)  and  gin)  are  so  called  quadrature  mirror  filters,  which  are  finite  impulse 
response  (FIR)  filters  that  allow  for  perfect  reconstmction  [3]  as  shown  in  Figure  2.7. 

Figure  2.7  illustrates  a  one- step  analysis  and  synthesis  process,  where  c^.^j  and  dJ_^^  are 
the  approximation  coefficients  and  detail  coefficients,  respectively.  The  lowpass  and 
highpass  filters  are  represented  by  h  and  g,  respectively,  and  i  2  represents  the 
decimation  by  two  operation.  Mallat’ s  lowpass  and  highpass  filters  are  simply  the 
scaling  basis  function  and  wavelet  basis  function  found  in  Equation  (2.10).  This 
decomposition  can  be  expanded  into  a  multi- scale  filter  bank  stmcture,  as  shown  in 
Figure  2.8,  where  N  refers  to  the  decomposition  level  and  “WL”  indicates  the  Wavelet 
analysis  operation.  Figure  2.8  demonstrates  the  approximation  coefficients  from  the 
previous  scale  are  used  to  form  the  detail  and  approximation  coefficients  of  the  following 
scale.  Note  that  the  number  of  scales  is  not  limited  to  three  as  shown  in  Figure  2.8,  but  is 
limited  by  the  length  of  the  data  segment.  The  result  is  that  the  frequency  spectmm  is 
sequentially  cut  in  half  as  shown  in  Figure  2.4  (a).  The  highest  resolution  coefficients 
labeled  “Detail  N  =  1”  in  Figure  2.8  represent  the  highpass  filtering  signal  contribution. 

Then,  the  next  resolution  takes  the  remaining  lower  half  of  the  spectmm  and  splits  it  in 

two.  By  using  the  factor  of  2,  Mallat  formed  what  are  commonly  referred  to  as  dyadic 

11 


scales.  Hence,  orthonormal  bandpass  filters  are  formed  via  MaUat’s  theorem  using  only 
lowpass  and  highpass  filtering  operations.  From  Equation  (2.1 1)  the  approximation 
coefficients  of  the  previous  scale  are  used  to  form  the  approximation  and  detail 
coefficients  of  the  next  higher  scale,  and  the  process  continues  until  the  desired  scale  is 
reached. 


Analysis  Synthesis 

Figure  2.7:  MaUat’s  filter  bank  wavelet  analysis  and  synthesis.  After  [3]. 


Figure  2.8:  Multi-level  wavelet  decomposition. 

Similar  to  the  STFT,  the  resulting  DWT  coefficients  can  be  represented  in  two- 
dimensional  time  and  frequency.  The  resulting  time- scale  visualization  is  called  the 
scalogram. 

Figure  2.9  presents  the  scalogram  of  the  chirp  from  Equation  (2.3),  to  be  visually 

compared  to  Figure  2. 1 .  Decomposition  level  one  shows  the  presence  of  high  frequency 
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during  the  second  half  of  the  signal.  The  resolution  of  high  frequency  components  was 
chosen  to  be  low  by  convention  since  it  is  assumed  that  the  majority  of  the  observed 
signals  reside  in  the  lower  half  of  the  spectrum.  For  the  case  where  most  of  the  signal 
hes  in  the  upper  half  of  the  spectmm,  a  different  multiresolution  technique  can  be  em¬ 
ployed  where  the  higher  frequencies  have  highest  resolution.  In  this  work  it  is  assumed 
that,  due  to  fast  processing  and  sampling  speeds,  most  observed  signals  will  he  in  the 
lower  half  of  the  spectmm.  Hence,  MaUat’s  multi-level  filter  decomposition  as  discussed 
above  is  applied  as  the  method  of  choice  throughout  this  work. 

The  sinusoidal  based  signal  such  as  the  chirp  in  Equation  (2.3)  is  more  easily  in¬ 
terpreted  as  a  constant  amphtude  linear  chirp  by  analyzing  the  spectrogram.  Hence,  the 
spectrogram  and  scalogram  are  often  used  in  concert  to  correctly  identify  the  characteris¬ 
tics  of  a  broad  class  of  signals. 


Time  (n) 

Figure  2.9:  Linear  chirp  Scalogram. 

The  previous  chapter  introduced  MaUat’s  theory  and  the  use  of  MRA  as  a  filter 
bank  stmcture.  The  filter  bank  spUts  the  signal  spectmm  into  dyadic  resolutions  with  the 
highest  resolution  representing  the  upper  half  of  the  spectmm.  When  a  signal  is  present. 
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signal  energy  is  found  in  one  or  more  of  the  frequency  scales.  Noise  energy  can  be 
present  throughout  each  of  the  dyadic  frequency  scales.  The  following  chapter  intro¬ 
duces  signal  processing  techniques  developed  to  clean  up  noisy  information. 
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m.  SIGNAL  COMPOSITION 


A  variety  of  modulation  formats  are  used  in  modem  communications,  which  leads 
to  a  large  variety  of  transmitted  and  received  signals.  Although  signals  come  in  various 
time- varying  shapes,  sizes,  frequencies,  and  phases,  signal  processing  techniques  have 
been  developed  to  digitally  format  signals  for  lossless  reconstmction. 


A.  SIGNALS 

Deterministic  functions  convey  information,  which  can  be  transmitted  in  the  form 
of  a  signal  using  the  electromagnetic  spectmm.  Analog  or  digital  signals  are  transmitted 
from  one  point  in  space  to  another  with  the  intention  of  being  received  and  analyzed  at 
the  receiver.  The  formation  of  a  signal  with  the  least  redundancy  is  based  on  the  func¬ 
tions’  required  bandwidth  as  defined  by  the  samphng  theorem  [6,9]. 

In  order  to  form  a  discrete  signal  from  a  continuous  analog  signal,  the  analog  sig¬ 
nal  is  sampled  at  the  samphng  interval  T  by  multiplying  by  an  impulse  train: 


(0  =  /(O  £  5  (t  -  nr). 

n=-oo 


(3.1) 


The  samphng  interval  T  is  chosen  smah  enough  so  the  sampled  data  can  be  used  to 
reconsfruct  the  original  signal.  Intuitively,  the  fewer  data  points  available  and  the  higher 
the  signal’s  frequency  the  more  difficult  it  wih  be  to  estimate  the  original  signal. 

Transforming  Equation  (3.1)  into  the  frequency  domain  by  noting  that  multiphca- 
hon  in  the  time  domain  is  equivalent  to  convolution  in  the  frequency  domain  leads  to  the 
foUowing  expression: 
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which  further  evaluates  as  fohows: 
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The  sampled  spectrum  found  in  Equation  (3.3)  contains  F((o)  with  periodic  identical 
rephcas  of  itself.  The  period,  which  is  related  to  the  sampling  rate  =  l/T ,  is  important 
because  if  it  is  too  small  then  overlap  and  aliasing  will  occur.  Figure  3.1  (a)  shows  the 
case  where  the  sampling  frequency  is  large  enough,  and  the  original  F((o)  can  be  fuUy 
recovered  with  a  lowpass  filter.  Figure  3.1  (b)  illustrates  an  undesirable  case  where  the 
original  signal  cannot  be  recovered  due  to  the  ahasing.  Notice  that  the  original  signal  can 
be  recovered  as  in  Figure  3.1  (a)  when  the  sampling  frequency  &  sampled  -  2^  1^  is  twice 
the  signal  bandwidth.  This  special  frequency  is  called  the  Nyquist  frequency  and  results 
in  the  minimum  safe  sampling  rate  for  real  signals  [11,12]. 


Figure  3.1: 
frequency; 


(0 


Sampled  frequency  spectmm:  (a)  Sampling  frequency  above  Nyquist 
(b)  Samphng  frequency  below  Nyquist  frequency. 


(0 


In  this  work  suitable  signals  are  constmcted  using  the  sampling  theorem  in  order 
to  test  various  denoising  algorithms.  The  following  section  describes  the  method  used  in 
this  work  to  form  a  sine  wave  with  frequency  set  to  a  uniform  random  variable  between  0 
and  71  radians. 

Equation  (2.3)  is  used  with  frequency  f^n)  as  a  uniform  random  variable  be¬ 
tween  0  and  4000  Hz.  Hence,  with  a  samphng  frequency  set  to  8000  Hz  we  meet  the 
requirements  of  the  Nyquist  rate.  The  input  signal  for  each  sinusoidal  simulation  was 
therefore  a  sine  wave  with  frequency  set  to  a  value  between  0  and  7i  radians  for  the 
duration  of  the  simulation.  Note  that  the  chirp  from  Equation  (2.3)  also  met  the  require¬ 
ments  of  the  Nyquist  Rate;  however  it  only  covered  frequencies  up  to  7t /2  radians. 
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B.  NOISE 


1.  Introduction 

Received  signals  contain  noise  over  the  entire  received  spectrum,  while  the  de¬ 
sired  signals  may  be  contained  in  a  small  dynamic  region  of  time  and  frequency.  Filters 
are  used  to  lower  the  signal- to- noise  ratio  (SNR)  by  eliminating  noise  from  the  parts  of 
the  spectmm  not  occupied  by  signal  energy.  Throughout  this  thesis  it  is  assumed  the 
signal  channel  has  Additive  White  Gaussian  Noise  (AWGN)  with  a  variance  equal  to 
one.  This  is  not  a  severe  restriction  since  pre- whitening  filters  can  be  employed  when  the 
noise  is  colored  [8,13]. 

2.  Noise  Power  Estimation 

When  noise  power  is  not  known  a  priori,  then  it  must  be  determined  using  noise 
power  estimation  techniques.  Once  the  noise  power  is  determined,  then  the  data  can  be 
normalized.  Two  techniques,  a  time -domain  technique,  and  a  wavelet- domain  technique 
are  discussed  in  this  section. 

The  time- domain  technique  rehes  on  determining  a  section  of  the  data  with  rela¬ 
tively  signal  free  characteristics.  This  section  of  data,  which  contains  predominately 
noise  energy  is  used  to  form  an  estimate  of  the  noise  power  for  the  data.  Once  an  esti¬ 
mate  of  the  noise  power  is  determined  it  is  normalized  to  unity,  by  dividing  the  entire 
signal  by  the  square  root  of  the  noise  power.  Then  the  resulting  signal  can  be  passed  as  a 
whole  or  in  smaller  subsections  through  the  proposed  denoising  algorithm. 

The  second  method  used  to  measure  the  noise  power  utilizes  knowledge  about  the 
wavelet  coefficients.  Mallat  [3:  p.  459]  shows  that  any  AWGN  at  the  input  to  the  wavelet 
transform  is  transformed  “at  the  finest  scale”  to  AWGN  of  the  same  variance.  If  it  can  be 
assumed  the  signal  lies  predominantly  in  the  lower  half  of  the  signal  specfrum,  then  the 
detail  coefficients  at  the  first  level  of  decomposition  will  be  primarily  AWGN.  Conse¬ 
quently,  measuring  the  variance  of  these  coefficients  and  normalizing  the  signal  as 
illustrated  above  also  leads  to  noise  normalized  data. 

In  either  case,  the  result  is  a  signal  combined  with  AWGN  with  a  variance  of  one. 


The  thesis  flow  diagram  of  Figure  1.1  shows  that  the  next  step  is  to  mn  the  signal  though 
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various  wavelet  and  time-based  denoising  processes.  The  following  chapter  introduces 
wavelet- based  denoising  concepts. 
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IV.  DENOISING  METHODS 


We  now  have  the  tools  necessary  to  delve  into  several  areas  of  denoising.  This 
work  focuses  primarily  on  (1)  thresholding  as  a  means  of  denoising  in  the  wavelet 
domain  due  to  some  important  characteristics  present  in  wavelet  coefficients;  and  (2) 
denoising  using  the  minimum  Mean  Squared  Error  (MSE)  technique  of  Wiener  EUtering 
in  the  time  domain.  We  also  introduce  a  hybrid  denoising  method.  Additionally,  trans¬ 
lation-invariant  denoising  is  explored  using  the  same  techniques  as  those  adopted  in  the 
orthogonal  wavelet  domain. 

A.  WAVELET  DOMAIN 

The  noisy  input  signal  can  be  thought  of  as  the  sum  of  the  desired  signal  compo¬ 
nent  (or  tme  signal)  and  the  Additive  White  Gaussian  Noise  (AWGN)  with  variance  of 
one: 

x(n)=  s(n)  -I-  n(n).  (4.1) 

desired  signal  noise  component 

It  has  been  shown  that  when  the  wavelet  basis  selected  is  well  matched  to  the  signal 
characteristics,  very  few  of  the  wavelet  detail  coefficients  are  influenced  by  the  signal, 
while  most  of  them  are  influenced  by  the  noise.  Therefore,  an  expression  for  the  wavelet 
coefficients  at  each  decomposition  level  can  be  described  by: 

yj(i)=  w  fi)  +  n.{i).  (4.2) 

desired  coefficients  noise  component 

In  addition,  the  desired  signal  coefficients  are  expected  to  be  of  larger  magnitude 
when  the  SNR  is  not  too  small.  Therefore,  denoising  is  accomplished  by  thresholding 
wavelet  coefficients,  thereby  eliminating  noise- only  coefficients  and  keeping  the  desired 
signal  coefficients  for  reconstmction. 

Upon  reconstmction,  the  MSE  normalized  by  the  signal  power  is  apphed  in  the 
following  manner  for  determining  how  closely  our  de- noised  signal  compares  to  the 
original  noiseless  signal: 
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(4.3) 


MSE  - 

normalized  y  \?-,  ’ 

E[s(n)  ] 

where  £[•]  is  the  expectation  operator,  and  s(n)  is  the  de-noised  signal.  A  visual 
comparison  can  be  made  between  x(n)  and  5(n)  when  the  original  signal  s(n)  is  not 
available.  The  rest  of  the  chapter  defines  and  illustrates  the  various  wavelet  thresholding 
techniques  currently  available  in  the  hterature. 

1.  Thresholding 

Thresholding  is  an  important  concept  in  denoising  and  compression,  because  as 
previously  stated,  a  few  detail  coefficients  hold  the  signal  information  when  the  wavelet 
basis  selected  is  well  matched  to  signal  characteristics,  while  the  effect  of  AWGN  on  the 
signal  is  the  same  over  all  the  coefficients  at  each  scale.  Note  that  the  approximation 
coefficients  that  do  not  contain  signal  energy  often  do  not  reside  at  or  near  zero,  as  do 
their  parent  detail  coefficients.  Hence,  thresholding  schemes  will  be  limited  to  detail 
coefficients. 

The  goal  of  this  section  is  to  define  threshold  levels  that  can  be  used  to  “kill”  de¬ 
tail  coefficients  at  each  decomposition  level  that  are  likely  to  contain  noise  energy  [1,2]. 
To  illustrate  the  above  properties  and  understand  their  usefulness,  it  is  helpful  to  compare 
and  contrast  the  coefficients  Wj  (/)  obtained  with  no  additive  noise  and  the  coefficients 

yj  (/)  obtained  with  noise.  Figures  4.1  (a)  and  (b)  provide  the  detail  coefficients  from  top 

to  bottom  starting  with  the  first  level  of  decomposition.  Figure  4.1  (a)  are  the  coefficients 
used  to  create  the  scalogram  in  Figure  2.9.  Note  the  similarities  between  the  two  figures. 
Figure  4.1  (a)  provides  a  plot  of  the  magnitude  of  the  desired  detail  coefficients  Wj  (i)  for 

a  linear  chirp,  and  Figure  4.1  (b)  shows  the  noisy  coefficients  (i) .  Notice  that  the  first 
level  of  decomposition  coefficients  in  Figure  4.1  (a)  are  zero  from  n  =  0  to  250,  hence 

the  thresholded  output  of  any  good  thresholding  scheme  will  be  close  to  zero  in  that 
region.  A  horizontal  fine  is  drawn  as  shown  in  Figure  4.1  (b)  and  labeled  T  for  threshold 
across  the  top  of  the  highest  coefficient  at  each  level  of  decomposition  produced  only  by 
noise.  The  height  of  the  fine  would  be  the  maximum  necessary  threshold  to  remove  all 
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the  noise  present  in  noise- only  contributed  coefficients.  Note  that  noise  is  not  removed 
from  signal  plus  noise  coefficients.  This  maximum  threshold  is  called  the  Universal  or 
Visual  threshold  and  will  be  described  in  more  detail  in  the  following  section. 


(a)  (b) 

Figure  4. 1 :  Magnitude  of  detail  coefficients  at  level  N  from  chirp  in  Equation  (2.3) 

with  signal  to  noise  ratio  of  6  dB.  (a)  Tme  noiseless  coefficients;  (b)  Noisy 
coefficients  with  Threshold  =  T. 


a.  Universal  Threshold 

The  visual  threshold  gets  its  name  from  the  relatively  noiseless  feature  of 
s(n)  because  it  guarantees  the  removal  of  all  noise- only  coefficients.  However,  since 
this  thresholding  method  “kills”  the  greatest  number  of  coefficients,  it  also  may  mistak¬ 
enly  eliminate  the  most  signal.  Therefore,  this  threshold  involves  the  greatest  risk  of 
loosing  signal- containing  coefficients.  (The  term  risk  is  used  synonymously  with  MSE  in 
some  texts.) 


The  probabihty  that  the  universal  threshold  is  greater  than  the  magnitude 
of  every  noise- only  coefficient  is  equal  to  one.  Recall  we  have  selected  the  noise  to  be 
AWGN  with  zero- mean  and  a  variance  one;  thus  the  wavelet  coefficients  distribution 
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function  is  Gaussian.  Note  that  only  a  small  percental  of  coefficients  wiU  have  large 
magnitudes.  Hence,  a  small  percentage  of  the  noise  may  be  mistaken  for  signal  coeffi¬ 
cients,  provided  the  interference  is  constiuctive,  if  the  threshold  is  not  sufficiently  high. 

Similarly,  a  coefficient  containing  signal  energy  may  be  mistaken  for  noise  when  de¬ 
structive  interference  is  encountered. 

As  the  sample  size  increases,  the  probabihty  or  likelihood  of  encountering 
larger  magnitude  coefficients  increases.  Hence,  a  larger  threshold  is  required  to  ensure 
the  elimination  of  noise-only  coefficients,  as  n  becomes  larger  so  that: 

Pr  (max|n^.(/)|  <rj  =  1,  (4.4) 

where  for  a  Gaussian  distribution  the  threshold  T  is  given  by: 

r=a^21n(n).  (4.5) 

Note  that  the  above  threshold  value  represents  the  upper  bound  needed. 

Since  the  threshold  value  does  not  level  off  unless  large  data  sizes  are  used,  as  shown  in 
Figure  4.2,  it  is  not  recommended  for  data  sets  involving  less  than  several  thousand 
samples.  Finally,  note  that  each  successive  level  of  detail  coefficients  contains  half  the 
number  of  samples  found  at  the  previous  level.  Therefore,  while  the  universal  threshold 
is  effective  for  the  first  level  of  decomposition,  there  may  not  be  enough  samples  in  the 
next  for  it  to  be  apphed.  This  property  limits  the  number  of  decomposition  levels  when 
solely  applying  the  universal  threshold. 
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Figure  4.2:  Illustration  of  universal  threshold  calculation  of  Threshold  T  as  a  function 

of  n . 

b.  SURE  Threshold 

The  previous  threshold  can  be  thought  of  as  the  threshold  that  produces 
the  maximum  risk,  whereas  the  Stein’s  Unbiased  Risk  Estimate  (SURE)  threshold  results 
in  the  minimum  risk.  These  two  thresholding  techniques  can  be  thought  of  as  the  upper 
and  lower  bounds  in  thresholding. 

Donoho  and  Johnstone  [1,2]  first  proposed  the  SURE  algorithm  as  a  result 
of  the  work  done  by  Stein  on  the  unbiased  risk  estimate  of  multivariate  normal  estimation 
problems.  They  determined  that  the  detail  coefficients  from  each  level  of  decomposition 
could  be  modeled  using  independent  multivariate  estimation  problems.  At  each  scale, 
through  statistical  analysis  of  the  coefficients,  a  threshold  is  adaptively  selected  based  on 
the  quantity  of  signal  present  in  the  coefficients.  By  applying  SURE  [2]  over  a  range  of 
thresholds,  the  minimum  risk  is  determined. 

The  SURE  algorithm  is  only  effective  for  SNR’s  greater  than  0  dB 
because  it  reties  heavily  on  statistical  information  obtained  from  the  desired  coefficients 
Wj  (i) .  Hence,  another  thresholding  technique  such  as  the  universal  threshold  should  be 

adopted  if  a  sparse  signal  condition  is  found.  The  Hybrid  threshold  simply  makes  a 
choice  between  either  the  universal  Threshold  or  the  SURE  threshold  based  on  which 
threshold  is  smaller. 
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c.  Bayesian  Thresholding 

Assuming  that  we  can  estimate  the  distribution  function  of  the  signal  coef¬ 
ficients,  a  relatively  new  thresholding  technique  can  be  apphed.  This  technique  of 
Bayesian  thresholding  is  not  used  in  this  work  because  no  a  priori  signal  information  is 
assumed.  However,  Bayesian  thresholding  has  been  shown  to  outperform  both  universal 
and  SURE  thresholding  in  image  denoising,  since  natural  images  often  can  be  modeled 
with  a  heavy  tailed  Gaussian  distribution  [15,16]. 


d.  Hard  and  Soft  Thresholding 

The  first  step  in  the  denoising  process  was  to  obtain  the  wavelet  transform 
of  the  signal  x(n)  using  a  suitable  basis  function.  Then,  a  threshold  is  obtained  using 
one  of  the  above  thresholding  methods.  Once  an  appropriate  threshold  is  determined  we 
must  decide  how  to  apply  it.  This  work  discusses  the  hard  and  soft  thresholding  tech¬ 
niques  as  they  pertain  to  a  given  threshold. 

Hard  thresholding  zeroes  out,  or  shrinks  [1,2],  the  coefficients  that  have 
magnitudes  below  the  threshold,  and  leaves  the  rest  of  the  coefficients  unchanged  [5]: 


dfi) 


\djii)\>T 

\d,ii)\<T. 


(4.6) 


Soft  thresholding  extends  hard  thresholding  by  shrinking  the  magnitude  of 
the  remaining  coefficients  by  T ,  producing  a  smooth  rather  than  abmpt  transition  to  zero 
[5,14]: 


dj(i) 


rsign[d^.(i)](l  dj(i)  \-T) 

1  0 


I  dj(i)  \>T 
otherwise. 


(4.7) 


The  smooth  transition  to  zero  results  in  noticeably  fewer  artifacts  upon  reconsfiuction, 
especially  when  deahng  with  image  denoising  [14].  Hence,  soft  thresholding  is  generally 
better  for  denoising  due  to  its  inherent  smoothing,  whereas  hard  thresholding  is  better 
suited  for  data  compression. 
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In  either  case,  perfect  reconstruction  is  not  possible  since  some  of  the  sig¬ 
nal  components  are  thrown  away  with  the  undesired  noise.  Furthermore,  any  threshold¬ 
ing  technique  other  than  the  universal  threshold  wiU  preserve  some  of  the  noise- only 
coefficients. 

A  drawback  to  thresholding  is  that  noise  affecting  the  remaining  signal 
coefficients  is  not  removed.  This  work  explores  algorithms  that  attempt  to  clean  the 
remaining  coefficients  prior  to  reconstmction.  Finally,  note  that  the  orthogonal  wavelet 
transform  is  not  translation- invariant,  i.e.,  the  wavelet  coefficients  change  when  the 
signal  is  transformed  using  different  time  shifts.  The  next  section  introduces  this  impor¬ 
tant  concept  as  an  extension  of  the  orthogonal  wavelet  transform  denoising. 

2.  Translation  Invariance 

Translation- invariant  denoising  is  a  term  coined  by  Donoho  and  Coifman  [17,18], 
which  illustrates  a  wavelet- based  denoising  method  that  attempts  to  remove  the  negative 
attributes  associated  with  individual  translations.  There  are  two  important  properties 
associated  with  orthonormal  wavelet  coefficients  that  lead  to  further  developments  of 
wavelet  denoising  theory.  First,  wavelet  coefficients  represent  a  projection  of  a  signal 
into  its  signal  subspace.  Second,  they  are  translation  variant  in  that,  for  any  right  or  left 
shifted  input  signal  x(n-x ),  one  can  expect  a  shghtly  different  projection  into  signal 
subspace.  It  is  the  ahgnment  or  translation  combined  with  the  effect  of  additive  noise 
that  prevents  the  orthonormal  wavelet  from  obtaining  the  hue  signal  subspace. 

Any  one  projection  will  contain  artifacts  or  spurious  oscillations  caused  by  the 
alignment  of  a  discontinuity  in  the  signal  and  the  wavelet.  This  is  called  the  pseudo- 
Gibbs  phenomena  [8,17].  This  phenomenon  depends  on  the  signal  and  basis  function  and 
is  locahzed  to  a  small  percentage  of  wavelet  coefficients.  An  advantage  of  denoising 
using  the  wavelet  transform  vs.  the  Fourier  transform  is  that  the  Gibbs  phenomenon  is 
distributed  over  all  the  Fourier  domain  coefficients  whereas  it  is  localized  to  a  few 
coefficients  in  the  wavelet  domain  [17].  Intuitively,  we  can  expect  an  average  of  resul¬ 
tant  signal  subspace  projections  to  be  a  closer  approximation  to  the  hue  projection. 
Additionally,  depending  on  alignment,  the  discontinuities  due  to  singularities  may  be 
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sparse  or  frequent,  hence  an  average  tends  to  reduce  the  magnitude  of  the  artifacts 
overall. 

Since  each  similar  projection  contains  redundant  information  this  could  be 
thought  of  as  a  form  of  diversity.  As  with  diversity,  for  the  same  bit  rate,  the  more 
diversity  we  have,  the  slower  the  information  flow.  Since  it  is  desired  to  optimize  the 
process,  an  efficient  systematic  process  called  cycle- spinning  was  introduced  by  Donoho 
and  Coifman  [17,18]. 

3.  Cycle-Spinning 

Denoising  with  the  addition  of  cycle- spinning  utUizes  the  orthonormal  wavelet 
transform,  consequently  no  new  transformation  is  required.  The  process  of  cycle- spin¬ 
ning  simply  involves  right  shifting  and  left  shifting  by  a  preset  amount  prior  to  the  next 
level  of  decomposition,  and  shifting  back  and  averaging  prior  to  the  next  level  of  recon- 
stmction.  Figure  4.3  (a)  illustrates  the  hue  output  coefficients  at  each  level  of  decompo¬ 
sition  of  a  constant  amphtude  hnear  chirp  as  a  result  of  cycle- spinning.  Notice  that  each 
decomposition  level  has  redundant  information  in  the  form  of  redundant  signal  subspace 
projections.  The  first  decomposition  level  consists  of  the  2048  coefficients  from  6145  to 
8192  and  contains  only  two  projections,  each  of  which  consist  of  1024  coefficients.  Note 
the  first  level  detail  coefficients  of  the  orthonormal  wavelet  transform  consist  of  1024 
coefficients.  The  projections  are  a  result  of  the  first  step  of  cycle -spinning,  which  sepa¬ 
rately  passes  right  and  left  shifted  versions  or  translations  of  the  noise  free  signal  through 
the  orthonormal  wavelet  transform.  This  process  results  in  two  upper  projections  or 
detail  projections  and  two  lower  projections  or  approximation  projections.  Through  a 
recursive  process,  each  lower  projection  results  in  two  upper  and  lower  projections  in  the 
next  level  of  decomposition.  The  decomposition  levels  proceed  from  right  to  left  in 
Figure  4.3  such  that  the  second  decomposition  level  from  4097  to  6144  consists  of  four 
detail  projections.  Note  that  in  Figures  4.3  (b)  through  (d)  the  approximation  coefficients 
from  1  to  2048  are  not  thresholded,  as  is  the  case  with  orthonormal  thresholding. 

Figure  4.3  (b)  shows  the  result  of  passing  a  noisy  signal  through  the  translation- 
invariant  wavelet  transform.  Recall  that  each  projection  is  a  different  approximation  of 
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the  tnie  signal  subspace.  Therefore,  the  denoising  process  is  applied  independently  to 
each  projection  prior  to  averaging  and  reconstruction.  Figures  4.3  (c)  and  (d) 
demonstrate  the  outcome  of  soft  thresholding  with  the  visual  threshold  and  hybrid 
threshold  respectfully.  Note  that  since  the  hybrid  thresholding  algorithm  found 
significant  evidence  of  signal  presence  it  defaulted  to  the  SURE  method  of  thresholding 
rather  than  to  the  visual  threshold. 
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Figure  4.3:  Translation- invariant  denoising  of  constant  amphtude  linear  chirp  with 

SNR  =  9  dB  and  N  =  3.  (a)  Tme  wavelet  coefficients;  (b)  wavelet  coefficients  in 
AWGN;  (c)  wavelet  coefficients  de- noised  using  soft  visual  thresholding;  (d)  wavelet 
Coefficients  de- noised  using  soft  SURE  thresholding. 


In  order  to  enhance  the  denoising  process  further,  we  rely  on  the  convergence 
property  of  subspace  projections  [18].  This  property  suggests  that  the  projection  into  a 
subspace  becomes  closer  or  converges  toward  the  tme  projection  each  time  it  is  recur¬ 
sively  passed  through  the  cycle -spinning  process.  This  concept  is  explored  further  in  the 
following  section  on  recursive  cycle- spinning  [18]. 
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4. 


Recursive  Cycle -Spinning 


Figure  4.4  illustrates  the  process  of  passing  the  output  of  the  translation- invariant 
denoising  algorithm  through  the  identical  denoising  process  i  times.  This  process  can  be 
performed  recursively,  hence  the  term  recursive  cycle -spinning 


i,(M) 


■Si-i(n) 


x(n) 


Figure  4.4:  Illustration  of  denoising  using  multiple  passes. 


The  results  of  recursive  cycle- spinning  demonstrate  that  multiple  passes  result  in 
a  projection,  which  converges  toward  the  tme  projection  [18].  The  threshold  and  thresh¬ 
olding  method  will  play  a  large  part  in  the  abihty  of  the  algorithm  to  converge  to  the  tme 
projection,  however  the  number  of  passes  at  which  the  recursive  cycle -spinning  con¬ 
verges  to  the  tme  subspace  depends  on  the  level  of  the  threshold  because  of  inherent 
losses  the  signal  undergoes  in  the  thresholding  process.  Absolute  convergence  to  the  tme 
subspace  is  not  possible  because  of  the  loss  of  smaller  signal  components  through  the 
course  of  thresholding.  For  a  properly  chosen  threshold,  each  recursive  iteration  retains 
relatively  large  signal  energy;  however  a  large  portion  of  noise  is  removed.  The  projec¬ 
tions  converge  toward  the  signal  subspace,  and  at  the  same  time  away  from  the  noise 
subspace  at  the  end  of  each  recursive  iteration.  A  point  will  be  reached  where  the  re¬ 
moval  of  noise  at  each  level  of  decomposition  is  no  longer  the  dominating  factor.  At  this 
point  we  are  close  in  the  mean  squared  sense  to  the  tme  projection  and  any  additional 
iterations  may  prove  counterproductive.  Hence,  making  the  threshold  smaller  serves  to 
increase  the  number  of  productive  recursive  iterations.  The  “optimum”  threshold  will 
remove  the  most  noise  with  the  lowest  information  loss.  Note  that  this  work  does  not 
attempt  to  find  the  optimum  threshold  or  the  number  of  iterations  necessary  to  reach  the 
minimum  MSB.  We  take  an  intuitive  approach  to  picking  a  suitable  threshold  and 
sufficient  number  of  iterations  to  illustrate  the  usefulness  of  the  technique. 
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Although  hard  thresholding  preserves  the  most  signal  energy  remaining  in  the  co¬ 
efficients  prior  to  reconstmction,  it  provides  less  than  a  smooth  final  result.  Therefore, 
we  choose  to  apply  ten  recursive  iterations  using  hard  thresholding  followed  by  an 
additional  iteration  using  soft  thresholding.  This  method  keeps  more  of  the  signal  energy 
from  iteration  to  iteration,  while  providing  a  smoother  result. 

Intuitively,  the  threshold  that  requires  the  least  fewest  number  of  iterations  is  the 
universal  threshold  since  it  is  the  largest  threshold.  Further,  since  the  universal  threshold 
is  guaranteed  to  remove  aU  the  noise  on  the  first  iteration,  any  further  iterations  may 
prove  counterproductive  due  to  the  added  risk  of  signal  loss  from  the  higher  threshold. 

Using  the  same  argument  as  above,  the  threshold  that  requires  the  largest  number 
of  iterations  is  the  SURE  threshold  because  it  is  the  smallest  threshold  and  removes  the 
least  noise  and  signal  alike.  Recall  that  upon  reaching  a  certain  point,  any  additional 
iteration  will  cause  signal  degradation  as  a  result  of  the  additional  loss  of  signal  energy. 
This  effect  worsens  with  higher  thresholds,  since  they  impose  a  higher  risk.  Therefore,  to 
prevent  overshooting  convergence,  SURE  thresholding  is  adopted  as  a  method  of  choice 
for  this  work. 

The  convergence  abihty  of  the  algorithm  is  pointed  out  with  Eigures  4.5  (a) 
through  (d)  by  comparing  two  possible  recursive  cycle- spinning  arrangements  to  the 
noiseless  and  noisy  projection.  One  would  expect  that  as  long  as  there  is  no  overshoot, 
the  coefficients  in  the  case  where  the  most  iterations  are  performed  will  more  closely 
represent  the  noiseless  coefficients  seen  in  Eigure  4.5  (a).  Eigure  4.5  (c),  which  repre¬ 
sents  denoising  using  two  iterations  of  hard  thresholding  followed  by  one  iteration  of  soft 
thresholding  with  the  SURE  threshold  appears  to  be  further  from  the  tme  projection  than 
our  recursive  results  in  Eigure  4.5  (d).  Notice  that  the  lowest  decomposition  level  detail 
coefficients  of  Eigure  4.5  (d)  are  most  similar  to  their  counter  parts  of  Eigure  4.5  (a). 

The  convergence  of  the  upper  level  coefficients  therefore  appears  to  lag  the  lower 
coefficients  in  the  convergence  process.  Hence  convergence  overshoot  may  occur  first  at 
the  lowest  level  of  decomposition  and  produce  signal  degradation  of  low  frequency  signal 
components.  We  leave  this  determination  for  future  study. 
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Figure  4.5 :  Constant  amplitude  linear  chirp  in  translation- invariant  wavelet  domain 

prior  to  last  reconstruction:  (a)  Tme  noiseless  coefficients;  (b)  Coefficients  in  AWGN 
with  SNR  =  9  dB;  (c)  Recursive  cycle -spinning  with  two  hard  iterations  followed  by  one 
soft  iteration  using  SURE  threshold;  (d)  Recursive  cycle- spinning  with  ten  hard  iterations 
followed  by  one  soft  iteration  using  SURE  threshold. 


This  method  will  be  compared  to  other  denoising  schemes  using  the  MSE.  Addi¬ 
tional  clean  up  of  the  remaining  coefficients  wiU  also  be  attempted  using  the  time  domain 
techniques  discussed  in  the  next  section. 


B.  TIME  AND  WAVELET  DOMAIN 

Time- domain  filtering  techniques  such  as  Wiener  Entering  have  been  adopted  to 
minimize  the  mean  squared  error  of  the  filtered  signal.  These  proven  time- domain 
techniques  wiU  be  compared  to  wavelet- domain  thresholding  schemes,  described  previ¬ 
ously  in  this  chapter.  In  addition  to  direct  comparison,  we  consider  a  combination  of 
both  wavelet  and  time  domain  techniques  to  further  improve  the  denoising  process. 

The  following  section  describes  Wiener  filtering  principles.  The  Wiener  filter  it¬ 
self  is  not  used  since  we  would  require  a  priori  knowledge  of  the  signal  or  an  accurate 
calculation  of  AWGN  variance.  Although  this  thesis  makes  the  assumption  that  the  noise 
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statistics  are  known  or  can  be  determined,  we  select  Wiener  prediction  as  a  useful  signal 
processing  alternative. 

1.  Wiener  Filtering 

Wiener  filtering  is  well-known  in  the  literature  [13].  When  the  observation  is  of 
the  form: 


x(n)=  s(n)+r\(n),  (4.8) 

where  the  signal  s(n)  and  noise  ri(n)  are  independent,  the  equations  for  the  optimal  filter 
are  of  the  form: 
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where  R^(l)  =  R,(l)  +  and  R^  (/)  and  Rf^(l)  are  the  autocorrelation  function  for  the 
signal  and  noise,  respectively.  These  are  the  Wiener- Hopf  equations. 


In  order  to  predict  the  signal  (i.e.,  estimate  s{n  -I- 1)  instead  of  s{n)) ,  the  Wiener- 
Hopf  equations  change  shghtiy: 
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a.  Elimination  of  Predictor  Edge  Distortion 

The  adopted  method  consists  of  flipping  the  data  from  front  to  back,  con¬ 
volving  with  the  same  prediction  filter  from  above,  and  then  flipping  back  again.  Hence, 
two  representations  of  the  data  signal  have  been  determined  where  one  is  forward  pre¬ 
dicted  and  the  other  reverse  predicted.  The  last  two  accurate  samples  from  the  reverse 
prediction  are  put  in  their  respective  place  in  the  forward  predicted  data  as  shown  in 
Figure  4.6  to  eliminate  the  edge  effect  cause  by  the  forward  predictor.  Additionally,  the 
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last  two  samples  of  the  forward  predicted  data  are  kept,  as  shown  in  Figure  4.6.  Note  that 
h  is  the  same  in  both  cases;  however  the  predicted  values  are  different.  Figure  4.6 
illustrates  the  process  demonstrating  that  the  remaining  samples  not  affected  by  the  edge 
distortion  are  then  averaged  together. 
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Figure  4.6:  Resulting  predicted  data  using  forward  and  reverse  prediction. 
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There  are  two  distinct  advantages  found  when  using  the  above  method. 
First,  the  size  of  the  output  data  will  be  the  same  as  that  of  the  input  data.  Next,  the  edges 
will  be  smooth  when  the  segments  are  put  back  together  because  the  end  points  are  fuUy 
predicted  with  each  of  the  filter  coefficients.  These  properties  are  advantageous  when  a 
large  stream  of  data  is  to  be  analyzed  in  small  incremental  segments.  Thus,  a  long  data 
set  can  be  spht  up  into  smaller  segments  and  when  reassembled  will  not  have  noticeable 
edges. 

For  a  non- stationary  signal  such  as  the  constant  amplitude  hnear  chirp,  a 
windowed  version  of  Wiener  Prediction  as  described  above  is  implemented  on  individual 
segments  of  the  data  stream.  If  the  segment  consists  of  zeroes,  then  the  predictor  returns 
aU  zeroes  to  prevent  the  formation  of  a  singular  matrix.  Note  that  the  wavelet- domain 
coefficients  contain  AWGN  as  previously  stated;  hence  Wiener  filtering  will  prove  useful 
in  denoising  the  remaining  coefficients. 

2.  Wavelet-  and  Time -Domain  Techniques 

For  this  work  the  above  windowed  predictor  was  used,  and  the  most  suitable  data 
input  size  to  the  predictor  was  determined  through  experimentation.  For  this  work  the 
prediction  window  is  limited  to  sixteen  samples  for  first  level  decomposition  wavelet 
denoising  and  thirty- two  samples  for  time-based  denoising.  For  the  signals  considered, 

we  have  found  through  experimentation  that  the  smaller  sample  interval  is  necessary  in 
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the  wavelet  domain  for  best  results  when  applied  only  to  the  first  level  detail  coefficients. 
Figures  4.7  (a)  and  (b)  show  the  first  700  points  of  the  noiseless  and  noisy  constant 
amphtude  linear  chirp,  and  Figure  4.7  (c)  demonstrates  the  denoising  capabihty  achieved 
by  implementing  the  predictor  on  the  constant  amphtude  hnear  chirp.  Notice  there  is  no 
added  distortion  at  the  edges  of  the  waveform  segments. 

The  foUowing  section  implements  a  low  computational  cost  filtering  method  for 
comparison  purposes.  The  Median  Filter  performs  well  for  non- stationary,  but  low- 
frequency  signals  because  of  its  inherent  lowpass  filtering  characteristics. 

3.  Median  Filtering 

A  one- dimensional  median  filter  is  chosen  because  of  its  simphcity  and  low  com¬ 
putational  cost.  However,  it  requires  the  user  to  have  some  a  priori  knowledge  of  the 
signal  characteristics  because  it  does  not  perform  weU  if  the  signal  has  high  frequency 
components.  The  constant  amphtude  hnear  chirp  and  random  sinusoid  generated  in  this 
work  occupy  the  high  frequency  spectmm;  so  denoising  results  are  expected  to  be  poor. 
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Figure  4.7:  Constant  amphtude  hnear  chirp  in  time  domain  (First  700  samples):  (a) 

Tme  noiseless  signal;  (b)  Signal  in  AWGN  with  SNR  =  9  dB;  (c)  Time  domain  Wiener 
prediction  denoising  with  window  size  =128;  (d)  3^*^  order  Median  filter  denoising. 


33 


The  median  filter  or  medfiltl  in  Matlab  [20]  estimates  s(n)  by  looking  at  s(n  - 1), 
s(n),  and  s(n  +  1),  and  chooses  as  output  s(n)  the  median  value  of  those  3  values.  For 
segments  this  method  also  returns  the  same  number  of  samples  it  receives.  Figure  4.7  (d) 
demonstrates  the  median  filtering  denoising  ability  on  a  constant  amplitude  linear  chirp. 
From  visual  inspection  it  appears  to  do  fairly  well  at  denoising  the  low  frequency  portion 
of  the  signal. 

This  chapter  identified  several  wavelet-  and  time-domain  denoising  schemes.  The 
following  chapter  compares  the  preceding  denoising  scheme.  First,  orthonormal  wavelet 
thresholding  techniques  are  compared  to  the  time- domain  techniques,  which  are  followed 
by  a  separate  comparison  of  translation- invariant  thresholding  techniques.  Next,  we 
compare  the  best  wavelet  thresholding  schemes  from  the  orthonormal  wavelet  and 
translation  invariant  wavelet  transforms.  Finally,  we  compare  results  obtained  with  the 
Wiener  prediction  and  Median  filter  schemes. 


34 


V.  RESULTS 


A.  INTRODUCTION 

The  previous  chapters  discussed  some  important  principles  and  theory  which  form 
the  background  for  time -frequency  based  denoising.  The  following  section  compares 
performances  obtained  for  each  scheme  considered  in  this  work.  The  signals  used  are 
robust  in  that  they  encompass  a  broad  range  of  the  spectmm  and  no  a  priori  signal 
information  is  exploited. 

Throughout  the  following,  the  result  of  Wiener  prediction  on  the  data  in  the  time 
domain  is  used  as  a  benchmark.  Hence,  the  other  denoising  techniques  will  be  compared 
against  the  benchmark  set  by  time -domain  Wiener  prediction.  In  all  the  following 
simulations  the  results  from  one  hundred  distinct  simulations  are  averaged  in  order  to 
produce  a  more  precise  statistical  comparison. 

I.  Orthonormal  Wavelet  Denoising 

The  first  simulations  compare  and  contrast  orthonormal  wavelet  domain  tech¬ 
niques.  Figures  5.1  through  5.3  illustrate  MSB  results  between  wavelet-based  methods 
and  time-based  methods  for  the  constant  amphtude  linear  chirp  of  Equation  (2.3),  and 
demonstrate  the  effect  of  adding  additional  levels  of  decomposition  N  where  N  =  2,  6, 
and  9,  respectively.  Based  on  MSE  performance,  the  time-domain  Wiener  prediction 
clearly  produced  the  best  results  for  this  signal  since  its  average  MSE  is  predominately 
less  than  each  of  the  other  methods.  Notice  that  for  A  =  6  and  9  the  MSE  for  the  Or¬ 
thonormal  soft  SURE  thresholding  performs  better  than  when  N  =  2.  This  result  is  to  be 
expected  since  more  frequency  ranges  are  being  cleaned.  Note  also  that  nine  levels  of 
decomposition  did  not  provide  noticeably  different  results  than  six  levels  of  decomposi¬ 
tion.  Hence,  a  point  of  diminishing  returns  is  reached  where  the  added  computations  (for 
implementing  further  levels  of  decomposition)  do  not  provide  added  benefit. 
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Recall  that  visual  thresholding  is  less  desirable  for  smaller  data  segments  because 
of  its  lower  performance  due  to  the  large  threshold  size  relative  to  the  number  of  data 
samples.  However,  the  benefits  of  Visual  thresholding  may  outweigh  those  of  SURE 
thresholding  for  longer  data  samples. 


Figure  5.1:  MSE  vs.  SNR  for  orthogonal  wavelet- based  denoising  of  constant  ampli¬ 

tude  linear  Chirp  from  Equation  (2.3)  with  N  =  2  (Average  of  100  simulations). 
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Figure  5.2:  MSB  vs.  SNR  for  orthogonal  wavelet- based  denoising  of  constant  ampli¬ 

tude  linear  Chirp  from  Equation  (2.3)  with  N  =  6  (Average  of  100  simulations). 

In  addition,  results  show  that  Median  filtering  produces  the  worst  MSE  results,  as 
expected,  since  the  Median  filter  acts  as  a  lowpass  filter.  Hence,  Median  filtering  does 
not  achieve  acceptable  results  for  the  nonstationary  multi- frequency  scale  signals  used. 
Median  filtering  may  produce  better  results  in  a  higher  sample  rate  environment  where 
changes  in  signal  statistics  are  relatively  insignificant  over  the  length  of  the  filter. 
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Figure  5.3:  MSB  vs.  SNR  for  orthogonal  wavelet- based  denoising  of  constant  ampli¬ 

tude  linear  Chirp  from  Equation  (2.3)  with  N  =  9  (Average  of  100  simulations). 

Figure  5.4  illustrates  one  trial  of  the  de- noised  signal  Sin)  and  the  original  signal 
sin)  to  demonstrate  the  visual  differences  found  in  the  various  denoising  schemes  for 
N  =  9.  These  results  are  consistent  with  those  found  in  Figure  5.3  and  illustrate  the 
difference  between  hard  and  soft  thresholding.  Results  show  better  MSB  performances 
are  obtained  with  the  SURE  soft  threshold  than  with  the  SURE  hard  threshold  for  this 
signal. 
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Figure  5.4:  First  700  samples  of  constant  amplitude  linear  Chirp  from  Equation  (2.3) 

using  orthogonal  wavelet- based  denoising  with  N  =  9  and  SNR  =  10  dB.  (a)  Time- 
domains  order  Wiener  Prediction;  (b)  Time- domain  Size  3  Median  Filter;  (c)  Soft 
visual  thresholding;  (d)  Hard  visual  thresholding;  (e)  Soft  SURE  thresholding;  (f) 

Hard  SURE  thresholding  (Original  noiseless  signal  is  shown  in  blue). 
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Figure  5.5  demonstrates  the  average  effectiveness  of  soft  SURE  wavelet  denois- 
ing  over  time  domain  Wiener  prediction  on  sinusoidal  signals.  Note  that  for  this  signal,  a 
decrease  in  the  MSE  is  illustrated  for  both  hard  and  soft  SURE  thresholding  compared  to 
the  Wiener  prediction  benchmark.  The  SURE  soft  threshold  MSE  was  shghtly  lower 
than  the  time  domain  Wiener  prediction  in  this  case,  whereas  for  the  constant  amphtude 
linear  chirp  it  was  shghtly  higher.  Thus,  the  overaU  MSE  performance  with  regard  to  the 
performance  of  wavelet  and  Wiener  prediction  on  both  signal  types  is  almost  identical. 

The  advantage  of  wavelet  thresholding  in  this  case  is  not  in  the  MSE  performance. 

Recall  that  the  best- suited  size  of  the  prediction  window  was  dependent  on  a  priori 
experimentation  whereas  wavelet  thresholding  did  not  require  any  a  priori  information 
other  than  segment  size  to  determine  the  number  of  levels  of  decomposition.  The  next 
Section  explores  the  effect  of  cycle -spinning  using  the  same  comparisons  as  above. 

2.  Translation-Invariant  Wavelet  Denoising 

In  this  Section  cycle -spinning  results  are  reported  and  compared  using  the  wavelet 
thresholding  techniques  discussed.  Eor  ease  of  comparison,  nine  decomposition  levels 
are  used  in  the  comparison.  Eigures  5.6  and  5.7  average  the  results  for  the  various  de¬ 
noising  schemes  on  the  constant  amphtude  chirp  and  random  sinusoid,  respectively. 

Results  for  the  constant  amphtude  chirp  are  similar  to  those  obtained  with  orthonormal 
wavelet  denoising;  however  a  few  differences  are  noted.  Eirst,  it  appears  that  translation 
invariance  as  apphed  produces  slightly  enhanced  denoising  capability  for  each  of  the 
thresholding  schemes.  In  addition,  hard  thresholding  with  the  visual  threshold  shows  a 
marked  improvement,  considering  the  smah  sample  space  involved. 
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Figure  5.5:  MSB  vs.  SNR  for  orthogonal  wavelet- based  denoising  of  sine  wave  with 

frequency  randomly  set  to  a  value  between  0  and  half  sampling  frequency  with  N  =  9 
(Average  of  100  simulations). 
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Figure  5.6:  MSB  vs.  SNR  for  translation- invariant  wavelet-based  denoising  of  con¬ 

stant  amplitude  linear  Chirp  from  Equation  (2.3)  with  N  =  9  (Average  of  100 
simulations). 
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Figure  5.7  that  shows  the  best  MSB  performances  are  obtained  with  the  soft 
SURE  threshold.  Next,  we  compare  the  soft  SURE  threshold  results  to  those  obtained 
with  cycle- spinning.  In  addition,  we  compare  combined  scheme  performances. 


Eigure  5.7:  MSB  vs.  SNR  for  translation- invariant  wavelet-based  denoising  of  sine 

wave  with  frequency  randomly  set  to  a  value  between  0  and  half  sampling  frequency  with 
N  =  9  (Average  of  100  simulations). 


3.  Combined  Time  and  Wavelet  Based  Denoising 

This  section  compares  the  best  orthonormal  and  translation- invariant  thresholding 
schemes,  and  provides  results  for  the  combined  schemes.  Recall  that  in  a  combined 
scheme  either  Wiener  prediction  or  median  filtering  is  apphed  to  the  wavelet  coefficients 
after  thresholding  and  prior  to  reconstmction.  This  technique  addresses  the  fact  that  no 
attempt  is  made  in  thresholding  to  eliminate  noise  in  the  remaining  coefficients.  Thresh¬ 
olding  continues  to  be  vital  in  the  denoising  process  because  it  is  able  to  completely  zero 
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out  a  portion  of  the  noise- only  eoeffieients  whereas  no  sueh  guarantee  exists  with  Wiener 
predietion  and  median  filtering.  ReeaU,  however,  that  with  Wiener  predietion  one  must 
make  eertain  the  autoeorrelation  matrix  is  not  singular,  whereby  ensuring  any  string  of 
eoeffieients  thresholded  to  zero  remain  zero.  Verifying  whether  the  matrix  is  iU-eondi- 
tioned  will  be  left  to  future  work;  however  in  sueh  a  ease  a  small  quantity  of  noise  may 
be  added  to  the  data  segment  prior  to  taking  its  autoeorrelation.  Additionally,  reeall  that 
for  this  work  the  time-based  teehniques  are  restrieted  to  the  detail  eoeffieients  from  the 
first  level  of  deeomposition. 

Figures  5.8  and  5.9  show  that  applying  Wiener  predietion  in  the  preseribed 
manner  does  provide  a  shght  advantage  over  denoising  utihzing  only  thresholding  for  the 
sine  wave  signal.  This  result  is  eneouraging  beeause  it  illustrates  that  additional 
denoising  ean  be  aeeomphshed  in  eonjunetion  with  thresholding.  Additionally,  it  shows 
that  translation- invariant  denoising  provides  added  denoising  eapability  over  that  of 
orthonormal  wavelet  denoising. 

Figure  5.10  provides  a  visual  representation  of  the  first  700  samples  of  the  de- 
noised  eonstant  amplitude  ehirp  s(n)  and  the  original  signal  s(n).  The  de- noised  signal 

in  Figure  5.10  (e)  demonstrates  the  effeetiveness  of  the  eombined  algorithm;  we  notiee 
the  relative  noise- free  nature  of  the  result.  Figure  5.10  (f)  shows  the  (huge)  degradation 
introdueed  by  applying  the  median  filter  to  the  eoeffieients  prior  to  reeonstmetion. 
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Mean  Squared  Error 


Figure  5.8:  MSB  vs.  SNR  for  combined  wavelet  and  time  denoising  of  constant 

amplitude  linear  Chirp  from  Equation  (2.3)  with  N  =  9  (Average  of  100  simulations). 
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Figure  5.9:  MSB  vs.  SNR  for  combined  wavelet  and  time  denoising  of  sine  wave  with 

frequency  randomly  set  to  a  value  between  0  and  half  sampling  frequency  with  N  =  9 
(Average  of  100  simulations). 
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Figure  5.10:  First  700  samples  of  constant  amplitude  linear  Chirp  from  Equation  (2.3) 
using  combined  wavelet  and  time-based  denoising  with  N  =  9  and  SNR  =  10  dB.  (a) 
Time  domain  3  order  Wiener  Prediction;  (b)  Time  domain  Size  3  Median  Filter;  (c) 
Orthonormal  soft  SURE  thresholding;  (d)  Translation- invariant  soft  SURE  thresholding; 
(e)  Translation- invariant  soft  SURE  thresholding  with  Wiener  Prediction  prior  to  recon¬ 
struction;  (f)  Translation- invariant  soft  SURE  thresholding  with  Median  filtering  prior 
to  reconstmction. 
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The  following  Section  provides  the  results  for  recursive  cycle -spinning  as  defined 
in  this  work  and  includes  a  variant  of  recursive  cycle -spinning  that  includes  Wiener 
prediction  just  prior  to  the  last  reconstruction  step. 

4.  Recursive  Translation-Invariant  Wavelet  Denoising 

This  section  combines  the  best  results  from  the  above  sections  and  compares  them 
to  recursive  cycle -spinning  and  a  combined  version  of  recursive  cycle- spinning  and 
Wiener  prediction. 

Figures  5.11  and  5.12  show  that  at  higher  SNR  levels  the  combined  recursive 
method  with  Wiener  prediction  performed  prior  to  the  last  reconstmction  has  the  best 
MSB  performance.  Notice  that  in  the  case  of  the  constant  amplitude  hnear  chirp  illus¬ 
trated  in  Figure  5.11,  the  combined  scheme  produces  better  results  at  lower  SNR’s  than 
for  the  random  sine  wave  illustrated  in  Figure  5.12.  Hence,  a  priori  information  about 
the  signal  is  necessary  in  the  implementation  of  this  combined  scheme;  therefore  it 
should  not  be  used  in  an  unknown  signal  environment. 

The  following  section  performs  the  denoising  algorithms  on  experimental  data 
sets.  No  a  priori  information  was  available  for  these  data  sets.  Due  to  the  large  sample 
rate  during  data  collection  these  segments  contain  a  large  number  of  samples. 
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Figure  5.11:  MSB  vs.  SNR  for  recursive  translation- invariant  wavelet  denoising  of 

constant  amplitude  linear  Chirp  from  Equation  (2.3)  with  N  =  9  (Average  of  100  simula¬ 
tions). 
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SNR(dB) 


Figure  5.12:  MSB  vs.  SNR  for  recursive  translation- invariant  wavelet  denoising  of  sine 
wave  with  frequency  randomly  set  to  a  value  between  0  and  half  sampling  frequency  with 
N  =  9  (Average  of  100  simulations). 
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5.  Denoising  of  Experimental  Data  Sets 

The  experimental  data  we  worked  with  are  referred  to  as  data  sets  A  through  F. 
They  are  de- noised  via  several  of  the  above  methods  as  shown  in  Figures  5.13  through 
5.26.  No  a  priori  information  has  been  collected  on  this  data,  hence  statistical  analysis 
was  performed  on  a  relatively  signal- free  portion  of  the  signal  to  determine  whether  it 
could  be  modeled  white  and  Gaussian.  Whiteness  of  the  noise  was  evaluated  by  com¬ 
puting  correlation  estimates  of  noise- only  data  segments,  while  the  Gaussian  assumption 
was  tested  with  quantile- quantile  plots.  Simulations  show  that  the  noise  segments  exhibit 
some  colored  properties  and  deviate  somewhat  from  the  Gaussian  assumptions.  The 
deviations  were  deemed  small  enough  so  that  we  could  stiU  attempt  to  consider  the  noise 
distortions  as  white  and  Gaussian  for  the  purpose  of  Wavelet  denoising.  An  additional 
pre- whitening  step  could  be  introduced  for  potentially  more  accurate  results  at  a  higher 
computational  cost.  Figures  5.13  through  19  compare  effects  of  the  visual  and  SURE 
threshold  and  the  orthogonal  transformation  to  the  translation- invariant  transformation. 

In  this  case  it  is  not  possible  to  notice  any  noticeable  difference  between  denois- 
ing  techniques  that  included  Prediction  and  those  that  did  not.  A  possible  explanation  is 
that  the  signal  did  not  change  significantly  during  the  prediction  window  duration  due  to 
the  large  samphng  rate.  Hence,  a  priori  information  about  the  signal  and  the  samphng 
frequency  is  necessary  to  correctly  apply  Wiener  prediction  in  both  time  and  wavelet 
domains. 

In  addition,  note  that  the  visual  thresholding  technique  essentially  eliminates  aU 
noise,  as  expected.  Since  the  data  segment  is  very  large  in  this  case  the  negative  effects 
experienced  with  the  smaller  data  sets  are  not  experienced  as  illustrated  in  Figures  5.13 
through  5.15.  Unfortunately,  visual  thresholding  stiU  finds  the  threshold  with  the  greatest 
risk  of  signal  loss.  For  instance,  the  latter  half  of  the  data  stream  in  Figure  5.19  may  have 
small  amphtude  signal  components  that  the  SURE  thresholding  technique  picks  up. 

Visual  thresholding  places  the  threshold  above  those  components  and  eliminates  them, 
however,  resulting  in  greater  jeopardy  of  thresholding  desired  signal  coefficients. 

Einally,  as  seen  by  comparing  Eigures  5.13  (e)  and  (f),  there  is  no  apparent 
difference  viewed  between  the  orthogonal  transformation  and  the  translation- invariant 
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transformation  for  these  data  sets.  Therefore,  the  added  cost  of  applying  the  translation- 
invariant  transform  does  not  justify  the  benefit  for  these  signals. 

Figures  5.20  through  26  illustrate  the  effect  of  recursive  cycle- spinning  on  the 
visual  and  SURE  soft  thresholding.  From  top  to  bottom  are  placed  the  original  data  set, 
Wiener  predicted  data  with  a  larger  window  size  of  512,  translation- invariant  soft  visual, 
recursive  cycle -spinning  with  one  iteration  of  hard  visual  followed  by  one  iteration  of 
soft  visual,  translation- invariant  soft  SURE,  and  recursive  cycle -spinning  with  ten 
iterations  of  hard  SURE  thresholding  followed  by  one  iteration  of  soft  SURE  threshold¬ 
ing. 

Figures  5.22  (e)  and  (f)  show  a  slight  denoising  advantage  of  recursive  cycle- spin¬ 
ning  over  its  translation- invariant  soft  SURE  counterpart.  In  most  cases  for  these  signals, 
however,  the  denoising  effect  of  recursive  cycle -spinning  was  unnoticeable  as  seen  in 
Figures  5.21  (e)  and  (f).  Thus,  the  cost  involved  in  recursive  cycle- spinning  is  much 
greater  than  its  benefit. 

This  chapter  systematically  identified  the  best  denoising  schemes  for  the  signals 
considered  and  applied  those  schemes  to  experimental  data.  The  following  chapter 
concludes  which  scheme  is  the  method  of  choice  for  the  signals  and  signal  lengths 
considered.  In  addition,  we  make  recommendations  for  further  study. 
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Figure  5.13:  Denoising  of  Data  segment  A.  (a)  Data  prior  to  denoising;  (b)  Time 
domain  3  order  Wiener  Prediction  with  window  size  of  32;  (c)  Orthonorml  soft  visual 
thresholding;  (d)  Translation- invariant  soft  visual  thresholding;  (e)  Orthonormal  soft 
SURE  thresholding;  (f)  Translation- invariant  soft  SURE  thresholding. 
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Figure  5.14:  Denoising  of  Data  segment  B.  (a)  Data  prior  to  denoising;  (b)  Time 
domain  3  order  Wiener  Prediction  with  window  size  of  32;  (c)  Orthonorml  soft  visual 
thresholding;  (d)  Translation- invariant  soft  visual  thresholding;  (e)  Orthonormal  soft 
SURE  thresholding;  (f)  Translation- invariant  soft  SURE  thresholding. 
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Figure  5.15:  Denoising  of  Data  segment  C.  (a)  Data  prior  to  denoising;  (b)  Time 
domain  3  order  Wiener  Prediction  with  window  size  of  32;  (c)  Orthonorml  soft  visual 
thresholding;  (d)  Translation- invariant  soft  visual  thresholding;  (e)  Orthonormal  soft 
SURE  thresholding;  (f)  Translation- invariant  soft  SURE  thresholding. 
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Figure  5.16:  Denoising  of  Data  segment  D.  (a)  Data  prior  to  denoising;  (b)  Time 
domain  3  order  Wiener  Prediction  with  window  size  of  32;  (c)  Orthonorml  soft  visual 
thresholding;  (d)  Translation- invariant  soft  visual  thresholding;  (e)  Orthonormal  soft 
SURE  thresholding;  (f)  Translation- invariant  soft  SURE  thresholding. 
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Figure  5.17:  Denoising  of  Data  segment  E.  (a)  Data  prior  to  denoising;  (b)  Time 
domain  3  order  Wiener  Prediction  with  window  size  of  32;  (c)  Orthonorml  soft  visual 
thresholding;  (d)  Translation- invariant  soft  visual  thresholding;  (e)  Orthonormal  soft 
SURE  thresholding;  (f)  Translation- invariant  soft  SURE  thresholding. 
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Figure  5.18:  Denoising  of  Data  segment  F.  (a)  Data  prior  to  denoising;  (b)  Time 
domain  3  order  Wiener  Prediction  with  window  size  of  32;  (c)  Orthonorml  soft  visual 
thresholding;  (d)  Translation- invariant  soft  visual  thresholding;  (e)  Orthonormal  soft 
SURE  thresholding;  (f)  Translation- invariant  soft  SURE  thresholding. 
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Figure  5.19:  Denoising  of  Data  segment  G.  (a)  Data  prior  to  denoising;  (b)  Time 
domain  3  order  Wiener  Prediction  with  window  size  of  32;  (c)  Orthonorml  soft  visual 
thresholding;  (d)  Translation- invariant  soft  visual  thresholding;  (e)  Orthonormal  soft 
SURE  thresholding;  (f)  Translation- invariant  soft  SURE  thresholding. 
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Figure  5.20:  Denoising  of  Data  segment  A.  (a)  Data  prior  to  denoising;  (b)  Time 
domain  3  order  Wiener  Prediction  with  window  size  of  512;  (c)  Translation- invariant 
soft  visual  thresholding;  (d)  Translation- invariant  recursive  with  one  hard  visual  fol¬ 
lowed  by  one  soft  visual ;  (e)  Translation- invariant  soft  SURE  thresholding;  (f)  Transla¬ 
tion-invariant  recursive  with  ten  hard  followed  by  one  soft  SURE  thresholding. 
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Figure  5.21 :  Denoising  of  Data  segment  B.  (a)  Data  prior  to  denoising;  (b)  Time 
domain  3  order  Wiener  Prediction  with  window  size  of  512;  (c)  Translation- invariant 
soft  visual  thresholding;  (d)  Translation- invariant  recursive  with  one  hard  visual  fol¬ 
lowed  by  one  soft  visual ;  (e)  Translation- invariant  soft  SURE  thresholding;  (f)  Transla¬ 
tion-invariant  recursive  with  ten  hard  followed  by  one  soft  SURE  thresholding. 
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Figure  5.22:  Denoising  of  Data  segment  C.  (a)  Data  prior  to  denoising;  (b)  Time 
domain  3  order  Wiener  Prediction  with  window  size  of  512;  (c)  Translation- invariant 
soft  visual  thresholding;  (d)  Translation- invariant  recursive  with  one  hard  visual  fol¬ 
lowed  by  one  soft  visual ;  (e)  Translation- invariant  soft  SURE  thresholding;  (f)  Transla¬ 
tion-invariant  recursive  with  ten  hard  followed  by  one  soft  SURE  thresholding. 
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Figure  5.23:  Denoising  of  Data  segment  D.  (a)  Data  prior  to  denoising;  (b)  Time 
domain  3  order  Wiener  Prediction  with  window  size  of  512;  (c)  Translation- invariant 
soft  visual  thresholding;  (d)  Translation- invariant  recursive  with  one  hard  visual  fol¬ 
lowed  by  one  soft  visual ;  (e)  Translation- invariant  soft  SURE  thresholding;  (f)  Transla¬ 
tion-invariant  recursive  with  ten  hard  followed  by  one  soft  SURE  thresholding. 
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Figure  5.24:  Denoising  of  Data  segment  E.  (a)  Data  prior  to  denoising;  (b)  Time 
domain  3  order  Wiener  Prediction  with  window  size  of  512;  (c)  Translation- invariant 
soft  visual  thresholding;  (d)  Translation- invariant  recursive  with  one  hard  visual  fol¬ 
lowed  by  one  soft  visual ;  (e)  Translation- invariant  soft  SURE  thresholding;  (f)  Transla¬ 
tion-invariant  recursive  with  ten  hard  followed  by  one  soft  SURE  thresholding. 
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Figure  5.25:  Denoising  of  Data  segment  F.  (a)  Data  prior  to  denoising;  (b)  Time 
domain  3  order  Wiener  Prediction  with  window  size  of  512;  (c)  Translation- invariant 
soft  visual  thresholding;  (d)  Translation- invariant  recursive  with  one  hard  visual  fol¬ 
lowed  by  one  soft  visual ;  (e)  Translation- invariant  soft  SURE  thresholding;  (f)  Transla¬ 
tion-invariant  recursive  with  ten  hard  followed  by  one  soft  SURE  thresholding. 
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Figure  5.26:  Denoising  of  Data  segment  G.  (a)  Data  prior  to  denoising;  (b)  Time 

domain  3  order  Wiener  Prediction  with  window  size  of  512;  (c)  Translation- invariant 
soft  visual  thresholding;  (d)  Translation- invariant  recursive  with  one  hard  visual  fol¬ 
lowed  by  one  soft  visual ;  (e)  Translation- invariant  soft  SURE  thresholding;  (f)  Transla¬ 
tion-invariant  recursive  with  ten  hard  followed  by  one  soft  SURE  thresholding. 
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VI.  CONCLUSIONS 


A.  METHOD  OF  CHOICE 

This  work  considered  and  compared  the  apphcation  of  time -based  and  wavelet- 
based  denoising  schemes  in  a  robust  signaling  environment.  A  time-based  predictor  was 
compared  against  time-based  median  filtering,  and  various  wavelet-based  techniques 
were  considered  along  with  a  combination  of  time-  and  wavelet-based  techniques.  A 
systematic  process  was  used  to  establish  which  would  be  the  method  of  choice  for  the 
signals  considered. 

First,  we  compared  time -based  techniques  with  the  orthonormal  wavelet  thresh¬ 
olding  techniques.  Results  show  that  the  soft  SURE  threshold  scheme  has  best  perform¬ 
ance  for  the  signal  types  and  length  considered.  Results  also  show  that  overall  the 
Wiener  prediction  scheme  produced  the  best  overall  results  for  the  two  signal  types 
considered.  Note,  however  that  Wiener  predictor  performance  is  window  size  dependent, 
and  Wiener  prediction  performance  became  similar  to  those  obtained  with  the  other 
schemes  for  an  appropriately  chosen  window  size.  In  addition,  the  results  show  that  the 
visual  thresholding  technique  was  inferior  to  the  SURE  thresholding  schemes  for  small 
data  sets.  However,  results  also  indicate  that  the  visual  threshold  performs  well  for  long 
data  sets  on  the  data  considered  in  this  study. 

Second,  we  compared  translation  invariance  soft  SURE  thresholding  schemes. 
Results  show  that  the  soft  SURE  thresholding  scheme  produced  the  best  results  for  the 
test  signals  considered.  Results  also  show  that  better  performances  are  obtained  for 
random  sinusoids  than  for  the  constant  amphtude  linear  chirp. 

Next,  we  compared  the  best  thresholding  schemes  obtained  for  the  orthonormal 
and  translation- invariant  wavelet  domain  against  the  combined  schemes.  We  found  that 
the  translation- invariant  soft  SURE  thresholding  outperformed  its  orthonormal  counter¬ 
part  consistently  over  the  range  of  SNR  levels  for  both  signal  sets.  In  addition,  Wiener 
prediction  apphed  to  the  thresholded  first  level  of  decomposition  detail  coefficients  just 
prior  to  reconstmction  added  shghtiy  to  the  denoising  abihty  for  SNR  (levels  from  two  to 
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eight  dB).  However,  results  show  no  noticeable  difference  for  our  test  data  for  the 
window  size  considered,  which  may  be  due  to  the  larger  data  sample  size. 

Finally,  we  compared  recursive  wavelet- based  techniques  to  the  best  performing 
orthonormal  and  translation- invariant  schemes.  We  found  that  the  recursive  scheme 
combined  with  Wiener  prediction  outperformed  other  schemes  for  SNR  levels  between 
six  and  eight  dB.  Thus,  results  show  that  the  signal  power  must  be  known  a  priori  in 
order  to  implement  such  a  technique.  In  addition,  the  Wiener  prediction  window  size 
must  be  optimized  for  the  sampling  environment. 

Our  results  suggest  that  the  translation- invariant  soft  SURE  thresholding  schemes 
should  be  employed  with  a  relatively  small  data  segment  and  no  a  priori  signal  informa¬ 
tion.  However,  for  larger  data  segments  soft  visual  thresholding  gives  the  greatest  noise- 
free  result,  yet  at  the  same  time  presents  the  greatest  risk  of  loosing  small  amphtude 
signal  components.  In  either  case  wavelet  thresholding  was  more  robust  than  time -based 
methods.  Therefore,  thresholding  is  indeed  robust  and  highly  capable  denoising  instm- 
ment.  Soft  SURE  thresholding  requires  the  least  quantity  of  a  priori  signal  information 
and  is  especially  robust  since  techniques  are  available  to  account  for  colored  noise 
environments. 

B.  RECOMMENDATION  FOR  FURTHER  STUDY 

Eurther  study  is  needed  in  the  area  of  recursive  cycle -spinning  and  Wiener  pre¬ 
dictive  denoising  of  coefficients.  Recall  that  Wiener  prediction  was  only  performed  on 
the  detail  coefficients  at  the  first  level  of  decomposition.  It  may  be  possible  to  implement 
an  adaptive  Wiener  prediction  scheme  that  chooses  its  window  size  based  on  statistics 
determined  from  a  range  of  window  sizes  for  each  level  of  decomposition.  In  addition,  it 
may  also  be  possible  to  de- noise  approximation  coefficients.  They  were  left  untouched  in 
our  study. 

The  visual  threshold  provides  the  largest  necessary  threshold  and  thereby  employs 
the  greatest  risk  of  signal  loss,  whereas  the  SURE  threshold  provides  the  minimum  risk 
of  signal  loss.  Recall  that  the  Hybrid  scheme  chooses  the  smaller  of  the  two  thresholds 
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where  the  SURE  threshold  is  based  on  evidence  of  signal  presence  and  the  visual  thresh¬ 
old  is  based  on  the  size  of  the  segment.  We  leave  for  further  work  another  version  of  the 
Hybrid  scheme  where  the  visual  threshold  is  chosen  if  it  is  less  than  the  SURE  threshold, 
and  an  average  of  the  two  thresholds  is  used  as  the  actual  threshold  when  the  SURE 
threshold  is  less  than  the  visual  threshold.  Hence  an  average  of  these  thresholds  may 
provide  a  good  compromise  between  the  two  schemes. 

Recursive  cycle- spinning  was  not  optimized  in  this  work,  however  the  results 
combined  with  Wiener  prediction  showed  promise.  Therefore,  optimization  of  the 
recursive  cycle- spinning  technique  for  each  level  of  decomposition  may  be  a  viable  area 
of  study. 
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APPENDIX  A  MATLABAVAVELAB  SOURCE  CODE 


The  Matlab  [20]  and  Wavelab  [10]  source  code  developed  and/or  modified  re¬ 
spectively  for  the  purposes  of  this  study  are  presented  in  this  appendix. 


A.  MAIN  PROGRAM  (AVERAGES  100  ITERATIONS) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

%  This  Matlab  Code  performs  averaging  of  100  separate  denoising 
%  iterations  for  N  wavelet  decomposition  levels . 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

%The  following  tines  of  code  are  used  to  create  the  constant 
%Linear  Chirp  from  Equation  (2.3).  Commented  out  if  Random  Sine 
%is  performed. 

% 

n=  1:2048; 

F=l:.9765:2000; 

Fs=8000; 

% 

%Create  Constant  amplitude  linear  chirp  from  0  to  one  fourth  sampling 
frequency 

% 

sig  1 =sin(2*pi.  *F./Fs.  *n) ; 
sig2=sigl/sqrt(mean(sigl  .^2)); 

% 

%Set  median  filter  size 

% 

medsize=3; 

% 

%Set  Wiener  Filter  Order 

% 

order=3; 

% 

%Set  number  of  decomposition  levels  to  be  performed  in  wavelet  routines 

% 

% 

forN=[2  6  9]; 

% 

for  nnn=l:100 
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% 

%Create  random  sinusoid  to  be  used  for  this  iteration  nnn 

% 

%Random  number  between  0  and  1 

F=rand; 

n=l:2048; 

Fs=2; 

% 

%Create  Random  Sinusoid  from  0  to  one  half  samphng  frequency 

% 

sig  1 =sin(2*pi*(F/Fs)  *n) ; 
sig2=sig  l/sqrt(mean(sig  1  .^2)) ; 

% 

%Store  snr  and  sigma  from  each  mn  in  Matrix  form 

% 

[L,snr(nnn,:),sigmal(nnn,:),sigma2(nnn,:),sigma3(nnn,:),... 
sigma4(nnn,:),sigma5(nnn,:),sigma6(nnn,:),orthohybridwcoef_soft]  =... 
OrthoDenoising(sig2, order, medsize,N); 

% 

end 

% 

figure(N); 

% 

subplot(  1,1, 1  ),semilogy  (mean(snr),mean(sigma  1  ),'bo- ') 

,hold  on,  semilogy  (mean(snr),mean(sigma2),'gx-') 

,hold  on,semilogy  (mean(snr),mean(sigma3),'r+-') 

,hold  on,semilogy  (mean(snr),mean(sigma4),'c*-') 

,hold  on,semilogy  (mean(snr),mean(sigma5),'ms-') 

,hold  on,semilogy  (mean(snr),mean(sigma6),'yd-') 

,hold  on,hold  off; 

% 

legend(sprintf('Time  Domain  %srd  order  Wiener  Prediction' ,num2str(order)  ),... 
spiintfCTime  Domain  Size  %s  Med  Filter',num2str(medsize)),... 

'Orthonormal  Soft  Visual  Thresholding','Orthonormal  Hard  Visual 
Thresholding',... 

'Orthonormal  Soft  SURE  Thresholding','Orthonormal  Hard  SURE  Thresholding'); 

% 

subplot(l,l,l),xlabel('SNR'),ylabel('Mean  Squared  Error'); 

% 

end 
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B.  DENOISING  FUNCTION  (ORTHONORMAL) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

%  This  Matlab  Function  calls  wavelab  and  modified  wavelab  denoising 
%  functions  [10].  Data  is  one  row  of  normalized  1  dimensional  data 
%  of  dyadic  length  2048=2^1 1 . 

%  Funetion  passes  the  signal  through  various 

%  Denoising  Sehemes  and  Passes  back  MSB  vs  SNR  for  each  of  those  schemes 

% 

%  Funetion  Syntax  as  follows: 

% 

%  [noisearray,sigmal,sigma2,sigma3,sigma4,sigma5,sigma6] 

%  =OrthoDenoising(sig2, order, medsize,N) 

% 

%  where: 

%  "sig2"  =  Tme  signal  passed  into  funetion 

%  "order"  =  Wiener  Predietor  Order 

%  "medsize"  =  median  filter  size 

%  "N"  =  Number  of  levels  of  decomposition  for  wavelet  transform 

%  "sigmal "  =  mean  squared  error  of  Time  Domain  Predietion 

%  "sigma2"  =  mean  squared  error  of  Time  Domain  Size  %s  Med  Filter 

%  "sigmaS"  =  mean  squared  error  of  Orthonormal  Soft  Visual  thresholding 

%  "sigma4"  =  mean  squared  error  of  Orthonormal  Hard  Visual  thresholding 

%  "sigmaS"  =  mean  squared  error  of  Orthonormal  Soft  SURE  thresholding 

%  "sigmab"  =  mean  squared  error  of  Orthonormal  Hard  SURE  thresholding 

%  "noisearray"  =  SNR  levels  [-10-6-3036  10]  dB 
% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

function  [noisearray,  sigmal,  sigma2,  sigma3,  sigma4,  sigmaS,  sigmab]... 
=OrthoDenoising(sig2,order,medsize,N) 

% 

%Make  a  eopy  of  signal 

% 

sigl=sig2; 

% 

%Eorm  array  of  ehosen  SNR's 

% 

noisearray=[-10  -6-3036  10]; 

% 

%oreate  noise  kernel 

% 

noise  =  randn(l,length(sigl)); 

% 

%find  noise  power 

% 
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noisepowe]^mean(noise.^2) ; 

% 

%normalize  noise  with  noisepower 

% 

noise  =  noise/sqrt(noisepower);%sets  noise  power  to  1 

% 

%For  loop  for  each  SNR;  Noise  seed  stays  the  same  for  each  SNR, 

%and  signal  power  changed  to  meet  requirements  of  noisearray. 

% 

for  count=l:length(noisearray) 

% 

%Change  the  signal  power  assuming  unit  variance  noise  for  desired  SNR 

% 

siglforsni^sqrt(  10^(.  1  *noisearray  (count)))  *sig  1 ; 

% 

%Add  normalized  noise  with  vai^l,  thereby  generating  SNR  dictated 

% 

si  =  noise  +  siglforsnr  ; 

% 

%Perform  Time  Domain  Prediction  Denoising 

% 

[FUteredPredict]  =  predict_window_time(sl, order); 

% 

%Perform  Time  Domain  Median  filtering 

% 

[FilteredMed]  =  medfiltl(sl,medsize); 

% 

%Perform  Orthonormal  Wavelet  Denoising 

% 

% 

%Create  quadrature  mirror  filter  for  Translation- invariant  Transform 

% 

QMF  =  MakeONFilter('Coiflef,5); 

% 

%Determine  length  of  signal  and  number  of  dyads  in  Data  using  Wavelab  routine 
%dyadlength  [10] 

% 

[lengthofsl,Jdyads]  =  dyadlength(sl); 

% 

%  Using  number  of  dyads  and  number  of  levels  of  decomposition 
%  find  degree  of  coarsest  scale,  ie.  if  N=6  and  Jdyads=5;  11-6=5 
% 

L=Jdyads-N; 

% 

%Perform  Orthonormal  Transformations  with  Shrinkage  using  modified 
%Wavelab  function  WaveShiink  [10]. 
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% 

[orthovisu_soft,  orthovisuwcoef_soft]  =  WaveShriiik_soft(sl,  'Visu',  L, 
QMF); 

[orthovisu_hard,  orthovisuwcoef_hard]  =  WaveShrirLk_hard(sl,  'Visu',  L, 
QMF); 

[orthosure_soft,  orthosurewcoef_soft]  =  WaveShriiik_soft(sl,  'SURE',  L, 
QMF); 

[orthosure_hard,  orthosurewcoef_hard]  =  WaveShriiik_hard(sl,  'SURE', 
L,  QME); 

% 

%Plot  results  of  the  various  filter  methods  over  the  top  of  the  noise  free  signal. 

% 

figure(N-l); 

subplot(6,l,l),plot(EilteredPredict(l:700),'r'),hold  on,  plot 
(siglforsnr(l:700),  'b'),  hold  off; 
subplot(6, 1 ,  l),ylabel('Amplitude') ; 
subplot(6,l,2),plot(PilteredMed(l:700),'r'),hold  on,  plot 
(siglforsnr(l:700),  'b'),  hold  off; 
subplot(6,l,2),ylabel('Amplitude'),title('(a)'); 
subplot(6,l,3),plot(orthovisu_soft(l:700),'r'),hold  on,  plot 
(siglforsnr(l:700),'b'),hold  off; 
subplot(6,l,3),ylabel('Amplitude'),title('(b)'); 
subplot(6,l,4),plot(orthovisu_hard(l:700),'r'),hold  on,  plot 
(siglforsnr(l:700),  'b'),hold  off; 
subplot(6,l,4),ylabel('Amplitude'),title('(c)'); 
subplot(6,l,5),plot(orthosure_soft(l:700),'r'),hold  on,  plot 
(siglforsnr(  1:700),  'b'),hold  off; 
subplot(6,l,5),ylabel('Amplitude')  ,title('(d)'); 
subplot(6,l,6),plot(orthosure_hard(l:700),'r'),hold  on,  plot 
(siglforsnr(l:700),'b'),hold  off; 
subplot(6, 1 ,6),ylabel('Amplitude'),  xlabel('(f)'),title('(e)'); 

% 

%  Perform  Mean  Squared  Error  Calculation. 

% 

%Time  Domain  Prediction 

[sigma  1  (count)]  =  msemed(PilteredPredict,siglforsnr); 

%Median  filtered 

[sigma2(count)]  =  msemed(PilteredMed,siglforsnr); 

%Orthonormal  Wavelet  Domain  Visual  Soft 
[sigma3(count)]  =  msemed(orthovisu_soft,siglforsnr); 

%  Orthonormal  Wavelet  Domain  Visual  Hard 
[sigma4(count)]  =  msemed(orthovisu_hard,siglforsnr); 

%  Orthonormal  Wavelet  Domain  SURE  Soft 
[sigma5(count)]  =  msemed(orthosure_soft,siglforsnr); 

%  Orthonormal  Wavelet  Domain  SURE  Hard 
[sigma6(count)]  =  msemed(orthosure_hard,siglforsnr); 
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end 


% 


C.  DENOISING  (TRANSLATION-INVARIANT) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

%  This  Matlab  Function  calls  wavelab  and  modified  wavelab  denoising 
%  functions  [10].  Data  is  one  row  of  normalized  1  dimensional  data 
%  of  dyadic  length  2048=2^1 1 . 

%  Function  passes  the  signal  through  various 

%  Denoising  Schemes  and  Passes  back  MSB  vs  SNR  for  each  of  those  schemes 

% 

%  Function  Syntax  as  follows: 

% 

%  [noisearray,sigmal,sigma2,sigma3,sigma4,sigma5,sigma6] 

%  =TransDenoising(sig2, order, medsize,N) 

% 

%  where: 

%  "sig2"  =  Tme  signal  passed  into  function 
%  "order"  =  Wiener  Predictor  Order 
%  "medsize"  =  median  filter  size 

%  "N"  =  Number  of  levels  of  decomposition  for  wavelet  transform 

%  "noisearray"  =  SNR  levels  [-10-6-3036  10]  dB 

%  "sigmal "  =  mean  squared  error  of  Time  Domain  Prediction 

%  "sigma2"  =  mean  squared  error  of  Time  Domain  Med  Filter 

%  "sigma3"  =  mean  squared  error  of  Translation- invariant  Soft  Visual  thresholding 

%  "sigma4"  =  mean  squared  error  of  Translation- invariant  Hard  Visual 

%  thresholding 

%  "sigmaS"  =  mean  squared  error  of  Translation- invariant  Soft  SURE  thresholding 

%  "sigma6"  =  mean  squared  error  of  Translation- invariant  Hard  SURE 

%  thresholding 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

function  [noisearray, sigmal, sigma2,sigma3,sigma4,sigma5,sigma6]  =... 
TransDenoising(sig2,order,medsize,N) 

% 

%Make  a  copy  of  signal 

% 

sigl=sig2; 

% 

%Eorm  array  of  chosen  SNR's 

% 

noisearray=[-10  -6-3036  10]; 

% 
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%create  noise  kernel 
% 

noise  =  randn(l,length(sigl)); 

% 

%fmd  noise  power 

% 

noisepowe]^mean(noise.^2); 

% 

%nomialize  noise  with  noisepower 

% 

noise  =  noise/sqrt(noisepower);%sets  noise  power  to  1 

% 

%For  loop  for  each  SNR;  Noise  seed  stays  the  same  for  each  SNR, 

%and  signal  power  changed  to  meet  requirements  of  noisearray. 

% 

for  count=l:length(noisearray) 

% 

%Change  the  signal  power  assuming  unit  variance  noise  for  desired  SNR 

% 

sig  1  forsni^sqrt(  10^( .  1  *noisearray  (count)))  *sig  1 ; 

% 

%Add  normalized  noise  with  vai^l,  thereby  generating  SNR  dictated 

% 

si  =  noise  +  siglforsnr  ; 

% 

%Perform  Time  Domain  Prediction  Denoising 

% 

[FilteredPredict]  =  predict_window_time(sl, order); 

% 

%Perform  Time  Domain  Median  filtering 

% 

[FilteredMed]  =  medfiltl(sl,medsize); 

% 

%Perform  Translation  Wavelet  Denoising 

% 

% 

%Create  quadrature  mirror  filter  for  Translation- invariant  Transform 

% 

QMF  =  MakeONFilter('Coiflet',5); 

% 

%Determine  length  of  signal  and  number  of  dyads  in  Data  using  Wavelab  routine 
%dyadlength  [10] 

% 

[lengthofsl,Jdyads]  =  dyadlength(sl); 

% 

%  Using  number  of  dyads  and  number  of  levels  of  decomposition 
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%  find  degree  of  coarsest  scale,  ie.  if  N=6  and  Jdyads=5;  1 1-6=5 
% 

L=Jdyads-N; 

% 

%Perform  Translation- invariant  Transformations  with  Shrinkage  using  modified 
%Wavelab  function  FWT_TI  [10]. 

% 

[transvisuwcoef_soft]=FWT_TI_visuthreshSoft(s  1  ,L,QMF); 
[transvisuwcoef_hard]=FWT_TI_visuthreshHard(s  1  ,L,QMF) ; 
[transsurewcoef_soft]=FWT_TI_SurethreshSoft(s  1  ,L,QMF); 
[transsurewcoef_hard]=FWT_TI_SurethreshHard(s  1  ,L,QMF); 

% 

%Perform  Inverse  Translation- invariant  Transform  with  Wavelab  function 
%IWT_TI  [10] 

% 

transvisu_soft=rWT_TI(transvisuwcoef_soft,QMF); 

transvisu_hard=rWT_TI(transvisuwcoef_hard,QMF); 

transsure_soft=rWT_TI(transsurewcoef_soft,QMF); 

transsure_hard=rWT_TI(transsurewcoef_hard,QMF); 

% 

%Plot  results  of  various  filter  methods  over  the  top  of  the  noise  free  signal. 

% 

figure(N-l); 

subplot(6,l,l),plot(FilteredPredict(l:700),'r'),hold  on,  plot  (siglforsnr(l:700), 
'b'),hold  off; 

subplot(6, 1 , 1  ),y  labelC  Amplitude') ; 

subplot(6,l,2),plot(FilteredMed(l:700),'r'),hold  on,  plot  (siglforsnr(l:700), 
'b'),holdoff; 

subplot(6,l,2),ylabel('Amplitude'),title('(a)'); 

subplot(6,l,3),plot(transvisu_soft(l:700),'r'),hold  on,  plot  (siglforsnr(  1:700), 
'b'),hold  off; 

subplot(6,l,3),ylabel('Amplitude'),titie('(b)'); 
subplot(6,l,4),plot(transvisu_hard(l:700),'r'),hold  on,  plot 
(siglforsnr(l:700),'b'),hold  off; 
subplot(6,l,4),ylabel('Amplitude'),title('(c)'); 

subplot(6,l,5),plot(transsure_soft(l:700),'r'),hold  on,  plot  (siglforsnr(l:700), 
'b'),holdoff; 

subplot(6, 1 ,5),ylabel(' Amplitude')  ,title('(d)'); 

subplot(6,l,6),plot(transsure_hard(l:700),'r'),hold  on,  plot  (siglforsnr(l:700), 
'b'),hold  off; 

subplot(6, 1 ,6),ylabel('Amplitude'),  xlabel('(f)'),titie('(e)'); 

% 

%  Perform  Mean  Squared  Error  Calculation. 

% 

%Time  Domain  Prediction 

[sigma  1  (count)]  =  msemed(FilteredPredict,siglforsnr); 
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%Translation- invariant  Wavelet  Domain  sure 
[sigma2(count)]  =  msemed(FilteredMed,siglforsnr); 

%Translation- invariant  Wavelet  Domain  sure 
[sigma3(count)]  =  msemed(transvisu_soft,siglforsnr); 

%Translation- invariant  Wavelet  Domain  Universal 
[sigma4(count)]  =  msemed(transvisu_hard,siglforsnr); 

%Translation- invariant  Wavelet  Domain  Universal 
[sigma5(count)]  =  msemed(transsure_soft,siglforsnr); 

%Time  domain  Med  Filtered 

[sigma6(count)]  =  msemed(transsure_hard,siglforsnr); 

% 

end 

D.  DENOISING  (COMBINED  TRANSLATION -INVARIANT  AND  WIENER 
PREDICTION) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

%  This  Matlab  Function  calls  wavelab  and  modified  wavelab  denoising 
%  functions  [10].  Data  is  one  row  of  normalized  1  dimensional  data 
%  of  dyadic  length  2048=2^11. 

%  Function  passes  the  signal  through  various 

%  Denoising  Schemes  and  Passes  back  MSB  vs  SNR  for  each  of  those  schemes 

% 

%  Function  Syntax  as  follows: 

% 

%  [noisearray  ,sigmal  ,sigma2,sigma3,sigma4,sigma5,sigma6] 

%  =CombinedDenoising(sig2, order, medsize,N) 

% 

%  where: 

%  "sig2"  =  Tme  signal  passed  into  function 
%  "order"  =  Wiener  Predictor  Order 

%  "medsize"  =  median  filter  size 

%  "N"  =  Number  of  levels  of  decomposition  for  wavelet  transform 

%  "noisearray"  =  SNR  levels  [-10-6-3036  10]  dB 

%  "sigmal "  =  mean  squared  error  of  Time  Domain  Prediction 

%  "sigma2"  =  mean  squared  error  of  Time  Domain  Med  Filter 

%  "sigma3"  =  mean  squared  error  of  Orthonormal  Wavelet  Domain  sure 
%  "sigma4"  =  mean  squared  error  of  Translation- invariant  Soft  SURE  thresholding 

%  "sigmaS"  =  mean  squared  error  of  Translation- invariant  Soft  SURE  thresholding 

%  with  Wiener  Prediction  on  1st  level  decomposition  detail  coeffs 

%  "sigma6"  =  mean  squared  error  of  Translation- invariant  Hard  SURE 
thresholding 

%  with  Median  Eiltering  on  1st  level  decomposition  detail  coeffs 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 
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function  [noisearray,sigmal,sigma2,sigma3,sigma4,sigma5,sigma6]  =... 
CombinedDenoising(sig2, order  ,medsize,N) 

% 

%Make  a  copy  of  signal 

% 

sigl=sig2; 

% 

%Form  array  of  chosen  SNR's 

% 

noisearray=[-10  -6-3036  10]; 

% 

%create  noise  kernel 

% 

noise  =  randn(l,length(sigl)); 

% 

%find  noise  power 

% 

noisepowe]^mean(noise.^2) ; 

% 

%normalize  noise  with  noisepower 

% 

noise  =  noise/sqrt(noisepower);%sets  noise  power  to  1 

% 

%For  loop  for  each  SNR;  Noise  seed  stays  the  same  for  each  SNR, 
%and  signal  power  changed  to  meet  requirements  of  noisearray. 

% 

for  count=l:length(noisearray) 

% 

%Change  the  signal  power  assuming  unit  variance  noise  for  desired  SNR 

% 

siglforsnr=sqrt(  10^(.  1  *noisearray(count)))*sigl ; 

% 

%Add  normalized  noise  with  vai^l,  thereby  generating  SNR  dictated 

% 

si  =  noise  -i-  siglforsnr  ; 

% 

%Perform  Time  Domain  Prediction  Denoising 

% 

[FilteredPredict]  =  predict_window_time(sl, order); 

% 

%Perform  Time  Domain  Median  filtering 

% 

[FilteredMed]  =  medfiltl(sl,medsize); 

% 

%Perform  Translation  Wavelet  Denoising 

% 
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% 

%Create  quadrature  mirror  filter  for  Translation- invariant  Transform 

% 

QMF  =  MakeONFilter('Coiflet',5); 

% 

%Determine  length  of  signal  and  number  of  dyads  in  Data  using  Wavelab  routine 
%dyadlength  [10] 

[lengthofsljdyads]  =  dyadlength(sl); 

% 

%  Using  number  of  dyads  and  number  of  levels  of  deeomposition 
%  find  degree  of  eoarsest  seale.  ie.  if  N=6  and  Jdyads=5;  1 1-6=5 

% 

L=Jdyads-N; 

% 

%Perform  Orthonormal, Translation- invariant,  and  eombined  Transformations 
%with  Shrinkage  using  modififieation  of  Wavelab  funetion  FWT_TI  [10], 

%and  our  predietion 
%function. 

% 

% 

%Perform  Orthonormal  Transformations  with  SURE  Shrinkage  using 
%modififieation  of 

%Wavelab  funetion  WaveShiink  [10]. 

% 

[orthosure_soft,  orthosurewcoef_soft]  =  WaveShrink_soft(sl,'SURE',L,QMF); 

% 

%Perform  Translation- invariant  Transformations  with  Shrinkage  using  modified 
%Wavelab  function  EWT_TI  [10]. 

% 

[transsurewcoef_soft]=EWT_TI_SurethreshSoft(s  1  ,E,QME); 
[transsurewcoef_soft_w]=EWT_TI_SurethreshSoft_weiner(sl,E,QME,order); 
[transsurewcoef_soft_m]=EWT_TI_SurethreshSoft_medfilt(s  1  ,E,QME) ; 

% 

%Perform  Inverse  Translation- invariant  Transform  with  Wavelab  function 
%IWT_TI  [10] 

% 

transsure_soft=rWT_TI(transsurewcoef_soft,QME); 

transsure_soft_w=rWT_TI(transsurewcoef_soft_w,QME); 

transsure_soft_m=rWT_TI(transsurewcoef_soft_m,QME); 

% 

%Plot  results  of  various  filter  methods  over  the  top  of  the  noise  free  signal. 

% 

figure(N-l); 

subplot(6,l,l),plot(EilteredPredict(l:700),'r'),hold  on,  plot  (siglforsnr(l:700),... 
'b'),hold  off; 

subplot(6, 1 , 1  ),y  labelC  Amplitude') ; 
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subplot(6,l,2),plot(FilteredMed(l:700),'r'),hold  on,  plot  (siglforsnr(l:700),... 

'b'),hold  off; 

subplot(6,l, 2), ylabelC  Amplitude'), title('(a)'); 

subplot(6,l,3),plot(orthosure_soft(l:700),'r'),hold  on,  plot  (sig lforsnr(  1:700),.. . 

'b'),hold  off; 

subplot(6,l,3),ylabel('Amplitude'),title('(b)'); 

subplot(6,l,4),plot(transsure_soft(l:700),'r'),hold  on,  plot  (siglforsnr(l:700),... 
'b'),holdoff; 

subplot(6,l,4),ylabel('Amplitude'),title('(c)'); 

subplot(6,l,5),plot(transsure_soft_w(l:700),'r'),hold  on,  plot  (siglforsnr(l:700),... 
'b'),hold  off; 

subplot(6,l,5),ylabel('Amplitude')  ,title('(d)'); 

subplot(6,l,6),plot(transsure_soft_m(l:700),'r'),hold  on,  plot  (siglforsnr(l:700),... 
'b'),hold  off; 

subplot(6, 1 ,6),ylabel('Amplitude'),  xlabel('(f)'),title('(e)'); 

% 

%  Perform  Mean  Squared  Error  Calculation. 

% 

%Time  Domain  Prediction 

[sigma  1  (count)]  =  msemed(FilteredPredict,siglforsnr); 

%Time  domain  Med  Filtered 

[sigma2(count)]  =  msemed(FilteredMed,siglforsnr); 

%Orthonormal  Wavelet  Domain  sure 
[sigma3(count)]  =  msemed(orthosure_soft,siglforsnr); 

%Translation- invariant  Soft  SURE  thresholding 
[sigma4(count)]  =  msemed(transsure_soft,siglforsnr); 

%Translation- invariant  Soft  SURE  thresholding  with  Wiener  Prediction 
[sigma5(count)]  =  msemed(transsure_soft_w,siglforsnr); 

%Translation- invariant  Soft  SURE  thresholding  with  Medfilter 
[sigma6(count)]  =  msemed(transsure_soft_m,siglforsnr); 

% 

end 

E.  DENOISING  (RECURSIVE,  SURE  SOFT  TRANSLATION -INVARIANT, 
AND  COMBINED) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

%  This  Matiab  Function  calls  wavelab  and  modified  wavelab  denoising 
%  functions  [10].  Data  is  one  row  of  normalized  1  dimensional  data 
%  of  dyadic  length  2048=2^1 1 . 

%  Function  passes  the  signal  through  various 

%  Denoising  Schemes  and  Passes  back  MSE  vs  SNR  for  each  of  those  schemes 

% 

%  Function  Syntax  as  follows: 

% 
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%  [noisearray,sigmal,sigma2,sigma3,sigma4,sigma5,sigma6] 

%  =Recursive(sig2, order, medsize,N) 

% 

%  where: 

%  "sig2"  =  True  signal  passed  into  function 
%  "order"  =  Wiener  Predictor  Order 

%  "medsize"  =  median  filter  size 

%  "N"  =  Number  of  levels  of  decomposition  for  wavelet  transform 

%  "noisearray"  =  SNR  levels  [-10-6-3036  10]  dB 

%  "sigmal"  =  mean  squared  error  of  Time  Domain  Prediction 

%  "sigma2"  =  mean  squared  error  of  Orthonormal  Wavelet  Domain  SURE  Soft 

%  "sigma3"  =  mean  squared  error  of  Translation- invariant  Soft  SURE  thresholding 

%  "sigma4"  =  mean  squared  error  of  Translation- invariant  Soft  SURE  thresholding 

%  with  Wiener  Prediction 

%  "sigmaS"  =  mean  squared  error  of  Recursive  cycle- spinning  10  hard  SURE 
%  Translation- invariant  followed  by  1  soft  SURE  Translation- invariant 

%  "sigma6"  =  mean  squared  error  of  Recursive  cycle -spinning  10  hard  SURE 

%  Translation- invariant  followed  by  1  soft  SURE  Translation- invariant 

%  with  Wiener  Prediction  on  1st  level  decomposition  detail  coeffs 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

function  [E, noisearray, sigmal, sigma2,sigma3,sigma4,sigma5,sigma6]  =... 
Recursive(sig2,order,medsize,N) 

% 

%Make  a  copy  of  signal 

% 

sigl=sig2; 

% 

%Eorm  array  of  chosen  SNR's 

% 

noisearray=[-10  -6-3036  10]; 

% 

%create  noise  kernel 

% 

noise  =  randn(l,length(sigl)); 

% 

%find  noise  power 

% 

noisepowe]^mean(noise.^2) ; 

% 

%normahze  noise  with  noisepower 

% 

noise  =  noise/sqrt(noisepower);%sets  noise  power  to  1 

% 

%Eor  loop  for  each  SNR;  Noise  seed  stays  the  same  for  each  SNR, 
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%and  signal  power  changed  to  meet  requirements  of  noisearray. 

% 

for  count=l:length(noisearray) 

% 

%Change  the  signal  power  assuming  unit  variance  noise  for  desired  SNR 

% 

sig  1  forsni^sqrt(  10^( .  1  *noisearray  (count)))  *sig  1 ; 

% 

%Add  normalized  noise  with  vai^l,  thereby  generating  SNR  dictated 

% 

si  =  noise  +  siglforsnr  ; 

% 

%Perform  Time  Domain  Prediction  Denoising 

% 

[FilteredPredict]  =  predict_window_time(sl, order); 

% 

%Perform  Time  Domain  Median  filtering 

% 

[FilteredMed]  =  medfiltl(sl,medsize); 

% 

%Perform  Translation  Wavelet  Denoising 

% 

% 

%Create  quadrature  mirror  filter  for  Translation- invariant  Transform 

% 

QMF  =  MakeONFilter('Coiflef,5); 

% 

%Determine  length  of  signal  and  number  of  dyads  in  Data  using  Wavelab  routine 
%dyadlength  [10] 

[lengthofsljdyads]  =  dyadlength(sl); 

% 

%  Using  number  of  dyads  and  number  of  levels  of  decomposition 
%  find  degree  of  coarsest  scale,  ie.  if  N=6  and  Jdyads=5;  1 1-6=5 

% 

L=Jdyads-N; 

% 

%Perform  Orthonormal, Translation- invariant,  and  combined  Transformations 
%with  Shrinkage  using  modifification  of  Wavelab  function  FWT_TI  [10],  and 
our  prediction 

%function. 

% 

% 

%Perform  Orthonormal  Transformations  with  Shrinkage  using  modified 
%Wavelab  function  WaveShiink  [10]. 

% 
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[orthosure_soft,  orthosurewcoef_soft]  = 

WaveShrink_soft(s  1  ,'SURE',L,QMF) ; 

% 

%Perform  Translation- invariant  Transformations  with  Shrinkage  using  modified 
%Wavelab  function  FWT_TI  [10]. 

% 

[transsurewcoef_soft]=FWT_TI_SurethreshSoft(s  1  ,L,QMF); 

[transsurewcoef_soft_w]=FWT_TI_SurethreshSoft_weiner(s  1  ,F,QMF, order) ; 

% 

%Perform  Inverse  Translation- invariant  Transform  with  Wavelab  function 
IWT_TI  [10] 

% 

transsure_soft=rWT_TI(transsurewcoef_soft,QMF); 

transsure_soft_w=rWT_TI(transsurewcoef_soft_w,QMF); 

% 

%Perform  Recursive  Cycle- spinning  with  Wavelab  function  rWT_TI  [10] 

% 

transsure_hard_r=s  1 ; 
for  n=l:10; 

% 

%Perform  Translation- invariant  Transformations  with  Shrinkage  using 

modified 

%  Wavelab  function  FWT_TI  [10]. 

% 

[transsurewcoef]=FWT_TI_Surethresh(transsure_hard_r,F,QMF); 

% 

%Perform  Inverse  Translation- invariant  Transform  with  Wavelab  function 

IWT_TI  [10] 

% 

transsure_hard_r=rWT_TI(transsurewcoef,QMF); 

% 

end 

% 

% 

%Perform  Translation- invariant  Transformations  with  Shrinkage  using  modified 
%Wavelab  function  FWT_TI  [10]. 

% 

[transsurewcoef]=FWT_TI_SurethreshSoft(transsure_hard_r,F,QMF); 

[transsurewcoef_wiener]=FWT_TI_SurethreshSoft_weiner(transsure_hard_r,F,Q 
MF, order); 

% 

%Perform  Inverse  Translation- invariant  Transform  with  Wavelab  function 
IWT_TI  [10] 

% 
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transsure_hard_r=rWT_TI(transsurewcoef,QMF); 

transsure_hard_wr=rWT_TI(transsurewcoef_wiener,QMF); 

% 

%Plot  results  of  various  filter  methods  over  the  top  of  the  noise  free  signal. 

% 

figure(N-l); 

subplot(6,l,l),plot(FilteredPredict(l:700),'r'),hold  on,  plot 
(siglforsnr(l:700),'b'),hold  off; 

subplot(6, 1 , 1  ),y  labelC  Amplitude') ; 
subplot(6,l,2),plot(orthosure_soft(l:700),'r'),hold  on,  plot 
(siglforsnr(l:700),'b'),hold  off; 

subplot(6, 1 ,2),ylabel('Amphtude'),title('(a)') ; 
subplot(6,l,3),plot(transsure_soft(l:700),'r'),hold  on,  plot 
(siglforsnr(l:700),'b'),hold  off; 

subplot(6,l,3),ylabel('Amplitude'),title('(b)'); 
subplot(6,l,4),plot(transsure_soft_w(l:700),'r'),hold  on,  plot 
(siglforsnr(l:700),'b'),hold  off; 

subplot(6,l,4),ylabel('Amplitude'),title('(c)'); 
subplot(6,l,5),plot(transsure_hard_r(l:700),'r'),hold  on,  plot 
(siglforsnr(l:700),'b'),hold  off; 

subplot(6,l,5),ylabel('Amplitude')  ,title('(d)'); 
subplot(6,l,6),plot(transsure_hard_wr(l:700),'r'),hold  on,  plot 
(siglforsnr(l:700),'b'),hold  off; 

subplot(6, 1 ,6),ylabel('Amplitude'),  xlabel('(f)'),title('(e)'); 

% 

%  Perform  Mean  Squared  Error  Calculation. 

% 

%Time  Domain  Prediction 

[sigma  1  (count)]  =  msemed(FilteredPredict,siglforsnr); 

%Orthonormal  Wavelet  Domain  SURE  Soft 
[sigma2(count)]  =  msemed(orthosure_soft,siglforsnr); 

%Translation- invariant  Soft  SURE  thresholding 
[sigma3(count)]  =  msemed(transsure_soft,siglforsnr); 

%Translation- invariant  Soft  SURE  thresholding  with  Wiener  Prediction 
[sigma4(count)]  =  msemed(transsure_soft_w,siglforsnr); 

%Recursive  cycle -spinning  10  hard  followed  by  1  soft 
[sigmaS (count)]  =  msemed(transsure_hard_r,siglforsnr); 

%Recursive  cycle -spinning  with  Wiener  Prediction 
[sigma6(count)]  =  msemed(transsure_hard_wr,siglforsnr); 

% 

end 

F.  MODIFIED  WAVELAB  FUNCTIONS 


I.  Orthonormal 
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a.  Modified  Waveshrink;  Waveshrink_hard.[10] 

%Modified  version  of  Waveshrink  from  Wavelab  [10]. 
function  [xh,  wcoefj  =  WaveShiink_hard(y,type,L,qmf) 

%  WaveShrink  --  Soft  Threshold  Shrinkage  Apphed  to  Wavelet 
%Coefficients 
%  Usage 

%  [xh,xwh]  =  WaveShrink(y,type,L,qnif) 

%  Inputs 

%  y  1-d  signal.  length(y)=  2^J 
%  Normalized  to  noise  level  1 !  (See  NoiseNorm) 

%  type  string.  Type  of  shrinkage  apphed: 

%  'Visu','SUREVHybrid','MinMaxVMAD' 

%  Optional;  default  ==  'Visu' 

%  L  Low-Frequency  cutoff  for  shrinkage  (e.g.  L=4) 

%  Should  have  L  «  J! 

%  qmf  Quadrature  Mirror  Filter  for  Wavelet  Transform 
%  Optional,  Default  =  Symmlet  8. 

%  Outputs 

%  xh  estimate,  obtained  by  applying  soft  thresholding  on 
%  wavelet  coefficients 
%  xwh  Wavelet  Transform  of  estimate 

% 

%  Description 

%  WaveShrink  smooths  noisy  data  presumed  to  have  noise  level  1 
%  by  transforming  it  into  the  wavelet  domain,  applying  soft 
%  thresholding  to  the  wavelet  coefficients  and  inverse  transforming. 

% 

%  The  theory  underlying  these  methods  is  described  in  a  variety  of 
%  papers  by  D.L.  Donoho  and  I.M.  Johnstone. 

% 

%  The  different  methods  of  selecting  thresholds  are  detailed 
%  in  their  articles. 

% 

%  See  Also 

%  FWT_PO,  rWT_PO,  MakeONFilter,  NoiseNorm,  RigorShrink 

% 

if  nargin  <  2, 
type  =  'Visu'; 
end 

if  nargin  <  3, 

L  =  3; 
end 

if  nargin  <  4, 

qmf  =  MakeONFilter('Symmlef,8); 
end 

% 
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[n,J]  =  dyadlength(y) ; 
wcoef  =  FWT_PO(y,L,qmf) ; 

% 

if  strcmp(type,'Visu'), 

wcoef((2^(L)+l):n)  =  VisuThresh(wcoef((2^(L)+l):n),’hard’) ; 
elseif  strcmp(type,'SURE'), 

wcoef  =  MultiSURE(wcoef,L); 
elseif  strcmp(type, 'Hybrid'), 

wcoef  =  MultiHybridhard(wcoef,L); 
elseif  strcmp(type,'MinMax'), 

wcoef((2^(L)+l):n)  =  MinMaxThresh(wcoef((2^(L)+l):n)) ; 
elseif  strcmp(type,'MAD'), 

wcoef  =  MultiMAD(wcoef,L); 
end 

% 

xh  =  IWT_PO(wcoef,L,qmf); 


% 

%  Copyright  (c)  1993-5.  Jonathan  Buckheit,  David  Donoho  and  Iain 
Johnstone 

% 


% 

%  Part  of  WaveLab  Version  802 

%  Built  Sunday,  October  3,  1999  8:52:27  AM 

%  This  is  Copyrighted  Material 

%  For  Copying  permissions  see  COPYING.m 

%  Comments?  e-mail  wavelab@stat.stanford.edu 

% 

b.  Modified  Waveshrink;  Waveshrink_soft[10]. 

function  [xh,  wcoef]  =  WaveShrink_soft(y,type,L,qmf) 

%  WaveShrink  --  Soft  Threshold  Shrinkage  Apphed  to  Wavelet 
%Coefficients 
%  Usage 

%  [xh,xwh]  =  WaveShrink(y,type,L,qmf) 

%  Inputs 

%  y  1-d  signal.  length(y)=  2^J 
%  Normalized  to  noise  level  1 !  (See  NoiseNorm) 

%  type  string.  Type  of  shrinkage  apphed: 

%  'Visu','SURE','Hybrid','MinMax','MAD' 

%  Optional;  default  ==  'Visu' 

%  L  Low-Frequency  cutoff  for  shrinkage  (e.g.  1^4) 

%  Should  have  L  «  J ! 

%  qmf  Quadrature  Mirror  Filter  for  Wavelet  Transform 
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%  Optional,  Default  =  Symmlet  8. 

%  Outputs 

%  xh  estimate,  obtained  by  applying  soft  thresholding  on 
%  wavelet  coefficients 
%  xwh  Wavelet  Transform  of  estimate 

% 

%  Description 

%  WaveShrink  smooths  noisy  data  presumed  to  have  noise  level  1 
%  by  transforming  it  into  the  wavelet  domain,  applying  soft 
%  thresholding  to  the  wavelet  coefficients  and  inverse  transforming. 

% 

%  The  theory  underlying  these  methods  is  described  in  a  variety  of 
%  papers  by  D.L.  Donoho  and  I.M.  Johnstone. 

% 

%  The  different  methods  of  selecting  thresholds  are  detailed 
%  in  their  articles. 

% 

%  See  Also 

%  FWT_PO,  rWT_PO,  MakeONFilter,  NoiseNorm,  RigorShrink 

% 

if  nargin  <  2, 
type  =  'Visu'; 
end 

if  nargin  <  3, 

L  =  3; 
end 

if  nargin  <  4, 

qmf  =  MakeONFilter('Symmlef,8); 
end 

% 

[n,J]  =  dyadlength(y) ; 
wcoef  =  FWT_PO(y,L,qmf) ; 

% 

if  strcmp(type,'Visu'), 

wcoef((2^(L)+l):n)  =  VisuThresh(wcoef((2^(L)+l):n)) ; 
elseif  strcmp(type,'SURE'), 

wcoef  =  MultiSUREsoft(wcoef,L); 
elseif  strcmp(type, 'Hybrid'), 

wcoef  =  MultiHybrid(wcoef,L); 
elseif  strcmp(type,'MinMax'), 

wcoef((2^(L)+l):n)  =  MinMaxThresh(wcoef((2^(L)+l):n)) ; 
elseif  strcmp(type,'MAD'), 

wcoef  =  MultiMAD(wcoef,L); 
end 

% 

xh  =  rWT_PO(wcoef,L,qmf); 
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% 

%  Copyright  (c)  1993-5.  Jonathan  Buckheit,  David  Donoho  and  Iain 
%Johnstone 

% 


% 

%  Part  of  WaveLab  Version  802 

%  Built  Sunday,  October  3,  1999  8:52:27  AM 

%  This  is  Copyrighted  Material 

%  For  Copying  permissions  see  COPYING.m 

%  Comments?  e-mail  wavelab@stat.stanford.edu 

% 


c.  Modified  MultSURE;  MultiSUREsoft  [10], 

function  ws  =  MultiSUREsoft(wc,L) 

%  MultiSURE  --  Apply  Shrinkage  to  Wavelet  Coefficients 
%  Usage 

%  ws  =  MultiSURE(wc,L) 

%  Inputs 

%  we  Wavelet  Transform  of  noisy  sequence  with  N(0,1)  noise 
%  L  low-frequency  cutoff  for  Wavelet  Transform 
%  Outputs 

%  ws  result  of  applying  SUREThresh  to  each  dyadic  block 

% 

[n,J]  =  dyadlength(wc); 
for  j=(J-l):-l:L 

wc(dyad(j))  =  SUREThreshsoft(wc(dyad(j))) ; 
end 

ws  =  we; 

% 

%  Copyright  (c)  1993-5.  Jonathan  Buckheit,  David  Donoho  and  Iain 
%Johnstone 

% 


% 

%  Part  of  WaveLab  Version  802 

%  Built  Sunday,  October  3,  1999  8:52:27  AM 

%  This  is  Copyrighted  Material 

%  For  Copying  permissions  see  COPYING.m 

%  Comments?  e-mail  wavelab@stat.stanford.edu 

% 
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d.  Modified  SUREThresh;  SUREThreshsoft  [10]. 
function  [x, thresh]  =  SUREThreshsoft(y) 

%  SUREThresh  --  Adaptive  Threshold  Selection  Using  Principle  of  SURE 
%  Usage 

%  thresh  =  SUREThresh(y) 

%  Inputs 

%  y  Noisy  Data  with  Std.  Deviation  =  1 
%  Outputs 

%  X  Estimate  of  mean  vector 
%  thresh  Threshold  used 

% 

%  Description 

%  SURE  referes  to  Stein's  Unbiased  Risk  Estimate. 

% 

thresh  =  ValSUREThresh(y); 

X  =  SoftThresh(y, thresh); 

% 

%  Copyright  (c)  1993-5.  Jonathan  Buckheit,  David  Donoho  and  Iain 
%Johnstone 

% 


% 

%  Part  of  WaveEab  Version  802 

%  Built  Sunday,  October  3,  1999  8:52:27  AM 

%  This  is  Copyrighted  Material 

%  Eor  Copying  permissions  see  COPYING.m 

%  Comments?  e-mail  wavelab@stat.stanford.edu 

% 

2.  Translation-Invariant 

a.  Modified  FWT_TI;  FWT_TI_visuthreshHard  [10]. 

function  wp  =  PWT_TI_visuthresh(x,E,qmf) 

%  PWT_TI  --  translation- invariant  forward  wavelet  transform 
%  Usage 

%  TIWT  =  PWT_TI(x,E,qmf) 

%  Inputs 

%  X  array  of  dyadic  length  n=2^J 
%  E  degree  of  coarsest  scale 
%  qmf  orthonormal  quadrature  mirror  filter 
%  Outputs 

%  TIWT  stationary  wavelet  transform  table 
%  formally  same  data  stiucture  as  packet  table 

% 
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%  See  Also 
%  IWT_TI 

% 


[n,J]  =  dyadlength(x); 

D  =  J-L; 

wp  =  zeros(n,D+l); 

X  =  ShapeAsRow(x); 

wp(:,l)  =  x'; 
ford=0:(D-l), 

forb=0:(2M-l), 
s  =  wp(packet(d,b,n),l)'; 
hsr  =  DownDyadHi(s,qmf); 
hsl  =  DowrLDyadHi(rshift(s),qmf); 

Isr  =  DowrLDyadLo(s,qmf); 

Isl  =  DowrLDyadLo(rshift(s),qmf); 
wp(packet(d+l,2*b  ,n),d+2)  =  VisuThresh(hsrVHard'); 
wp(packet(d+l,2*b+l,n),d+2)  =  VisuThresh(hsr,'Hard'); 
wp(packet(d+l,2*b  ,n),l  )  =  (Isr'); 
wp(packet(d+l,2*b+l,n),l  )  =  (Isl'); 
end 
end 


% 

%  Copyright  (c)  1994.  David  L.  Donoho 

% 


% 

%  Part  of  WaveLab  Version  802 

%  Built  Sunday,  October  3,  1999  8:52:27  AM 

%  This  is  Copyrighted  Material 

%  For  Copying  permissions  see  COPYING.m 

%  Comments?  e-mail  wavelab@stat.stanford.edu 

% 


b.  Modified  FWT_TI;  FWT_TI_visuthreshSoft  [10]. 

function  wp  =  FWT_TI_visuthresh(x,L,qmf) 

%  FWT_TI  --  translation- invariant  forward  wavelet  transform 
%  Usage 

%  TIWT  =  FWT_TI(x,L,qmf) 

%  Inputs 

%  X  array  of  dyadic  length  n=2^J 
%  L  degree  of  coarsest  scale 
%  qmf  orthonormal  quadrature  mirror  filter 
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%  Outputs 

%  TIWT  stationary  wavelet  transform  table 
%  formally  same  data  stmcture  as  packet  table 

% 

%  See  Also 
%  IWT_TI 

% 


[n,J]  =  dyadlength(x); 

D  =  J-L; 

wp  =  zeros(n,D+l); 

X  =  ShapeAsRow(x); 

wp(:,l)  =  x'; 
ford=0:(D-l), 

forb=0:(2M-l), 
s  =  wp(packet(d,b,n),l)'; 
hsr  =  DownDyadHi(s,qmf); 
hsl  =  DownDyadHi(rshift(s),qmf); 

Isr  =  DownDyadLo(s,qmf); 

Isl  =  DownDyadLo(rshift(s),qmf); 
wp(packet(d+l,2*b  ,n),d+2)  =  VisuThresh(hsrVSoff); 
wp(packet(d+l,2*b+l,n),d+2)  =  Visuthresh(hslVSoff); 
wp(packet(d+l,2*b  ,n),l  )  =  (Isr'); 
wp(packet(d+l,2*b+l,n),l  )  =  (Isf); 
end 
end 


% 

%  Copyright  (c)  1994.  David  L.  Donoho 

% 


% 

%  Part  of  WaveLab  Version  802 

%  Built  Sunday,  October  3,  1999  8:52:27  AM 

%  This  is  Copyrighted  Material 

%  For  Copying  permissions  see  COPYING.m 

%  Comments?  e-mail  wavelab@stat.stanford.edu 

% 


c.  Modified  FWT_TI;  FWT_TI_SUREthreshSoft  [10]. 

function  wp  =  FWT_TI_SurethreshSoft(x,L,qrnf) 

%  FWT_TI  --  translation- invariant  forward  wavelet  transform 
%  Usage 

%  TIWT  =  FWT_TI(x,L,qmf) 

93 


%  Inputs 

%  X  array  of  dyadic  length  n=2^J 
%  L  degree  of  coarsest  scale 
%  qmf  orthonormal  quadrature  mirror  filter 
%  Outputs 

%  TIWT  stationary  wavelet  transform  table 
%  formally  same  data  stmcture  as  packet  table 

% 

%  See  Also 
%  IWT_TI 

% 


[n,J]  =  dyadlength(x); 

D  =  J-L; 

wp  =  zeros(n,D+l); 

X  =  ShapeAsRow(x); 

% 

wp(:,l)  =  x'; 
ford=0:(D-l), 

forb=0:(2M-l), 
s  =  wp(packet(d,b,n),l)'; 
hsr  =  DownDyadHi(s,qmf); 
hsl  =  DownDyadHi(rshift(s),qmf); 
Isr  =  DownDyadLo(s,qmf); 

Isl  =  DownDyadLo(rshift(s),qmf); 
wp(packet(d+l,2*b  ,n),d+2)  = 
SUREThreshTransSoft(hsr') ; 
wp(packet(d+l,2*b+l,n),d+2)  = 
SUREThreshTransSoft(hsr) ; 
wp(packet(d+l,2*b  ,n),l  )  =  lsr'; 
wp(packet(d+  l,2*b+l,n),l  )  =  lsr; 
end 
end 


% 

%  Copyright  (c)  1994.  David  L.  Donoho 

% 


% 

%  Part  of  WaveLab  Version  802 

%  Built  Sunday,  October  3,  1999  8:52:27  AM 

%  This  is  Copyrighted  Material 

%  For  Copying  permissions  see  COPYING.m 

%  Comments?  e-mail  wavelab@stat.stanford.edu 

% 
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d.  Modified  FWT_TI;  FWT_TI_SUREthreshHard  [10]. 

function  wp  =  FWT_TI_SurethreshHard(x,L,qmf) 

%  FWT_TI  --  translation- invariant  forward  wavelet  transform 
%  Usage 

%  TIWT  =  FWT_TI(x,L,qmf) 

%  Inputs 

%  X  array  of  dyadic  length  n=2^J 
%  L  degree  of  coarsest  scale 
%  qmf  orthonormal  quadrature  mirror  filter 
%  Outputs 

%  TIWT  stationary  wavelet  transform  table 
%  formally  same  data  stmcture  as  packet  table 

% 

%  See  Also 
%  IWT_TI 

% 


[n,J]  =  dyadlength(x); 

D  =  J-L; 

wp  =  zeros(n,D-i-l); 

X  =  ShapeAsRow(x); 

wp(:,l)  =  x'; 
ford=0:(D-l), 

forb=0:(2M-l), 
s  =  wp(packet(d,b,n),l)'; 
hsr  =  DownDyadHi(s,qmf); 
hsl  =  DownDyadHi(rshift(s),qmf); 

Isr  =  DownDyadLo(s,qmf); 

Isl  =  DownDyadLo(rshift(s),qmf); 
wp(packet(d-i-l,2*b  ,n),d+2)  =  SUREThreshTrans(hsr'); 
wp(packet(d-i-l,2*b-i-l,n),d+2)  =  SUREThreshTrans(hsr); 
wp(packet(d-i-l,2*b  ,n),l  )  =  lsr'; 
wp(packet(d-i- l,2*b-i-l,n),l  )  =  lsr; 
end 
end 


% 

%  Copyright  (c)  1994.  David  L.  Donoho 

% 


% 

%  Part  of  WaveLab  Version  802 
%  Built  Sunday,  October  3,  1999  8:52:27  AM 
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%  This  is  Copyrighted  Material 
%  For  Copying  permissions  see  COPYING.m 
%  Comments?  e-mail  wavelab@stat.stanford.edu 

% 


e.  Modified  FWT_TI;  SUREThreshTrans  [10]. 
function  [x, thresh]  =  SUREThreshTrans(y) 

%  SUREThresh  --  Adaptive  Threshold  Selection  Using  Principle  of  SURE 
%  Usage 

%  thresh  =  SUREThresh(y) 

%  Inputs 

%  y  Noisy  Data  with  Std.  Deviation  =  1 
%  Outputs 

%  X  Estimate  of  mean  vector 
%  thresh  Threshold  used 

% 

%  Description 

%  SURE  referes  to  Stein's  Unbiased  Risk  Estimate. 

% 

thresh  =  ValSUREThresh(transpose(y)); 

X  =  HardThresh(y, thresh); 


% 

%  Copyright  (c)  1993-5.  Jonathan  Buckheit,  David  Donoho  and  Iain 
%Johnstone 

% 


% 

%  Part  of  WaveEab  Version  802 

%  Built  Sunday,  October  3,  1999  8:52:27  AM 

%  This  is  Copyrighted  Material 

%  Eor  Copying  permissions  see  COPYING.m 

%  Comments?  e-mail  wavelab@stat.stanford.edu 

% 


f  Modified  FWT_T1;  SUREThreshTransSoft  [10]. 

function  [x, thresh]  =  SUREThreshTransSoft(y) 

%  SUREThresh  --  Adaptive  Threshold  Selection  Using  Principle  of  SURE 
%  Usage 

%  thresh  =  SUREThresh(y) 

%  Inputs 
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%  y  Noisy  Data  with  Std.  Deviation  =  1 
%  Outputs 

%  X  Estimate  of  mean  vector 
%  thresh  Threshold  used 

% 

%  Description 

%  SURE  referes  to  Stein's  Unbiased  Risk  Estimate. 

% 

thresh  =  ValSUREThresh(transpose(y)); 

X  =  SoftThresh(y, thresh); 


% 

%  Copyright  (c)  1993-5.  Jonathan  Buckheit,  David  Donoho  and  Iain 
%Johnstone 

% 


% 

%  Part  of  WaveEab  Version  802 

%  Built  Sunday,  October  3,  1999  8:52:27  AM 

%  This  is  Copyrighted  Material 

%  Eor  Copying  permissions  see  COPYING.m 

%  Comments?  e-mail  wavelab@stat.stanford.edu 

% 

Combined 

a.  Modified  FWT_TI;  FWT_TI_SurethreshSoft_wiener  [10]. 

function  wp  =  PWT_TI_SurethreshSoft_wiener(x,E,qmf, order) 

%  PWT_TI  --  translation- invariant  forward  wavelet  transform 
%  Usage 

%  TIWT  =  PWT_TI(x,E,qmf) 

%  Inputs 

%  X  array  of  dyadic  length  n=2^J 
%  E  degree  of  coarsest  scale 
%  qmf  orthonormal  quadrature  mirror  filter 
%  Outputs 

%  TIWT  stationary  wavelet  transform  table 
%  formally  same  data  stmcture  as  packet  table 

% 

%  See  Also 
%  IWT_TI 

% 


[n,J]  =  dyadlength(x); 
D  =  J-E; 

wp  =  zeros(n,D-i-l); 
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X  =  ShapeAsRow(x); 

% 

wp(:,l)  =  x'; 
ford=0:(D-l), 

forb=0:(2M-l), 
s  =  wp(packet(d,b,n),l)'; 
hsr  =  DownDyadHi(s,qmf); 
hsl  =  DowrLDyadHi(rshift(s),qmf); 
Isr  =  DowrLDyadLo(s,qmf); 

Isl  =  DowrLDyadLo(rshift(s),qmf); 
size(SUREThreshTransSoft(hsr')) ; 
if  length(hsr')>length(x)/4 

wp(packet(d+l,2*b  ,n),d+2)  = 
predict_window(SUREThreshTransSoft(hsr'), order,  16); 
end 

if  length(hsr')<=length(x)/4 

wp(packet(d+l,2*b  ,n),d+2)  = 

SUREThreshTransSoft(hsr') ; 

end 

if  length(hsr)>length(x)/4 

wp(packet(d+l,2*b+l,n),d+2)  = 
predict_window(SUREThreshTransSoft(hsr), order,  16); 
end 

if  length(hsr)<=length(x)/4 

wp(packet(d+l,2*b+l,n),d+2)  = 

SUREThreshTransSoft(hsr) ; 

end 

wp(packet(d+l,2*b  ,n),l  )  =  lsr'; 
wp(packet(d+  l,2*b+l,n),l  )  =  lsr; 
end 
end 

% 

%  Copyright  (c)  1994.  David  L.  Donoho 

% 


% 

%  Part  of  WaveLab  Version  802 

%  Built  Sunday,  October  3,  1999  8:52:27  AM 

%  This  is  Copyrighted  Material 

%  For  Copying  permissions  see  COPYING.m 

%  Comments?  e-mail  wavelab@stat.stanford.edu 

% 
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b.  Modified  FWT_TI;  FWT_TI_SurethreshSoft_medfilt  [10]. 

limction  wp  =  FWT_TI_SurethreshSoft_medfilt(x,L,qmf, order) 

%  FWT_TI  --  translation- invariant  forward  wavelet  transform 
%  Usage 

%  TIWT  =  FWT_TI(x,L,qmf) 

%  Inputs 

%  X  array  of  dyadic  length  n=2^J 
%  L  degree  of  coarsest  scale 
%  qmf  orthonormal  quadrature  mirror  filter 
%  Outputs 

%  TIWT  stationary  wavelet  transform  table 
%  formally  same  data  stmcture  as  packet  table 

% 

%  See  Also 
%  IWT_TI 

% 

[n,J]  =  dyadlength(x); 

D  =  J-L; 

wp  =  zeros(n,D-i-l); 

X  =  Shape  AsRow(x); 

% 

wp(:,l)  =  x'; 
ford=0:(D-l), 

forb=0:(2M-l), 
s  =  wp(packet(d,b,n),l)'; 
hsr  =  DownDyadHi(s,qmf); 
hsl  =  DownDyadHi(rshift(s),qmf); 

Isr  =  DownDyadLo(s,qmf); 

Isl  =  DownDyadLo(rshift(s),qmf); 
size(SUREThreshTransSoft(hsr')) 

wp(packet(d-i-l,2*b  ,n),d+2)  = 
medfiltl(SUREThreshTransSoft(hsr'),3); 

size(wp(packet(d-i-l,2*b  ,n),d+2)) 

wp(packet(d-i-l,2*b-i-l,n),d+2)  = 
medfilt  1  (SUREThreshTransSoft(hsr),3) ; 

size(wp(packet(d-i- 1 ,2*b-i-l  ,n),d+2)) 

wp(packet(d-i-l,2*b  ,n),l  )  =  lsr'; 
wp(packet(d-i- l,2*b-i-l,n),l  )  =  lsr; 
end 
end 


% 

%  Copyright  (c)  1994.  David  L.  Donoho 
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% 


% 

%  Part  of  WaveLab  Version  802 

%  Built  Sunday,  October  3,  1999  8:52:27  AM 

%  This  is  Copyrighted  Material 

%  For  Copying  permissions  see  COPYING.m 

%  Comments?  e-mail  wavelab@stat.stanford.edu 

% 


G.  WINDOWED  PREDICTOR 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

%  This  Matlab  Function  performs  Windowed  Wiener  prediction 
%  to  the  given  "order"  on  sequencial  segments  of  Data. 

%  Input  must  be  2^Integer,  and  greater  than  W  in  length. 

%  Wiener  Prediction  conducted  forward  in  time  and  reverse  in  time 
%  The  Data  is  forward  predicted,  then  reverse 

%  predicted.  The  first  (order- 1)  data  bits  not  fuUy  covered  by  the  forward  filter 
%  are  replaced  by  the  data  bits  found  using  the  reverse  filter  and  vice  versa 
%  The  estimate  for  predicting  one  sample  ahead  is  s(n-i-l). 

% 

%  Function  Syntax  as  follows: 

% 

%  [Filtered]  =  predict_window(Data,order,W) 

% 

%  where: 

%  "Data"  =  True  signal  passed  into  function  length  2^Integer 

%  "order"  =  Wiener  Predictor  Order  or  predictor  size 

%  "W"  =  Window  Size 

%  "Filtered"  =  Predictor  output 

% 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

function  [Filtered]  =predict_window(Data,order,W) 

% 

%  Put  Data  in  proper  [1  length(Data)]  format 
Predictsize  =  size(Data); 
if  Predictsize(  1  )>  1 

Data=transpose(Data) ; 
end 

% 

%  Spht  Data  into  segnents  and  perform  denoising  on  each  segment 
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% 

for  g=0:(length(Data)AV)- 1 

% 

%Form  Rx  of  Wiener  Hopf  Equation 

% 

DataW=Data(l+g*W:W+g*W); 
Datasize=size(Data) ; 

DataW  size=size(DataW); 
AutoCorrX=xcorr(DataW) ; 

% 

%Form  Rx[0]..Rx[order- 1] 

% 

rx=AutoCorrX(W:W+order-  l)';%r(0),...r(order- 1) 

% 

%Fomi  Rs[l]..Rs  [order] 

% 

rs=AutoCorrX((W  + 1 ) :  W+order)'; 

% 

%Fmd  toeplitz  Rx 

% 

R=toeplitz(rx); 

% 

%Check  to  ensure  matrix  not  singular 

% 

flag=0; 

for  b=l:length(DataW) 

ifDataW(b)  ~=0 
flag  =1; 
end 
end 

if  flag  =1 

Rinv=inv(R); 

end 

% 

%If  aU  values  are  one  then  singular 
%Future  work  should  also  correct  for 
%ill- conditioned 

% 

if  flag  =  0 

Rinv=R; 

end 

% 

%Foim  prediction  Alter 

% 

h=Rinv*rs; 

% 
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%Convolve  S(n+1)  =  h[0]x[n]  +  h[l]x[n-l]  +  ...  +  h[p]x[n-p] 

% 

for  n=  1  :(W -  (order)) 

FF=DataW  (n:  (n+order- 1 )).  *fliplr(transpose(h)) ; 
Filtered(n+g*W)=sum(FF) ; 
end 

% 

%Ensure  data  still  in  correct  format 

% 

Filteredsize=size(Filtered) ; 
if  Filteredsize(  1  )>  1 

Filtered=transpose(Filtered) ; 
end 

% 

%Keep  the  filtered  values  that  touched  entire  filter 
%Hence  Keep  all  but  first  (order- 1)  values 

% 

Filteredl(l-i-g*W:W-i-g*W)=[DataW(l:order),  Filtered(l-i-g*W:W-i-  g*W- 
order )]; 

Filtered  1  size=size(Filtered  1 ) ; 

% 

%Ensure  data  stiU  in  correct  format 

% 

if  Filteredlsize(l)>l 

Filtered  1  =transpose(Filtered  1 ) ; 
end 

% 

%reverse  Data 

% 

DataW  l=fliplr(flipud(DataW)); 

% 

%Filter  in  reverse  direction 

% 

for  n=  1  :(W -  (order)) 

FFl=DataWl(n:(n-i-order-  l)).*fliplr(transpose(h)); 

Filtered3(n-i-g*  W)=sum(FF  1 ) ; 
end 

% 

%Ensure  data  stiU  in  correct  format 

% 

Filtered3size=size(Filtered3); 
if  FilteredS  size(  1  )>  1 

Filtered3=transpose(Filtered3); 

end 

% 

%Keep  the  filtered  values  that  touched  entire  filter 
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% 

Filtered2(l+g*W:W+g*W)=[DataWl(l:order),. . . 

Fi]tered3(  l+g*W  :W+g*W-order)] ; 

Filtered2size=size(Filtered2); 

% 

%Ensure  data  still  in  correct  format 

% 

if  Filtered2size(  1  )>  1 

Filtered2=transpose(Filtered2) ; 
end 

% 

%Flip  reversed  data  back,  such  that  it  compares  to  forward  data  point 
%by  point 

% 

Filtered2(  1  +g*W :  W+g*W)=fliplr(fhpud(Filtered2(  l+g*W :  W+g*W))); 

% 

%The  first  (order- 1)  data  bits  not  fuUy  covered  by  the  forward  filter 
%are  replaced  by  the  data  bits  found  using  the  reverse  filter  and  vice  versa 

% 

Filtered  1(1  -i-g*  W :  1  -i-g*  W  -i-order)=Filtered2(  1  -i-g*  W :  1  -i-g*  W  -i-order) ; 
Filtered2(W- (order-  l):W)=Filteredl(W- (order-  1):W); 
Filtered(l-i-g*W:W-i-g*W)=  .5*(Filtered2(l-i-g*W:W-i-  g*W)-i- 
Filteredl(l-i-g*W:W-i-g*W)); 

% 

%Ensure  data  stiU  in  correct  format 

% 

if  Predictsize(  1  )>  1 

Eiltered=transpose(Eiltered) ; 
end 

% 
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